NPDF5
NPDF5
develop reliability models for failure-time prediction under small failure-time sam-
model extends the work of Whitmore et al. 1998, to incorporate two new data-
time samples within dynamic environments where failure mechanisms are unknown,
there is a need for models that make use of auxiliary reliability information. In this
thesis we present models suitable for reliability data, where degradation variables
are latent and can be tracked by related observable variables we call markers.
and predictive inference equations for a data-structure that includes terminal ob-
more et. al. 1998 and show improvement in inference under small sample sizes. We
Bayesian support vector machine and discuss its place in degradation modeling. We
model a recurrent event process and failure-times. We compute the expected time to
failure using counting process theory and investigate the effect of the event process
by
Vasilis A. Sotiris
Dissertation Committee:
Professor Michael G. Pecht, Chair/Advisor
Professor Eric Slud, Chair/Co-Advisor
Professor Konstantina Trivisa
Professor Abhijit Dasgupta
Professor Peter Sandborn, Dean’s representative
UMI Number: 3495754
In the unlikely event that the author did not send a complete manuscript
and there are missing pages, these will be noted. Also, if material had to be removed,
a note will indicate the deletion.
UMI 3495754
Copyright 2012 by ProQuest LLC.
All rights reserved. This edition of the work is protected against
unauthorized copying under Title 17, United States Code.
ProQuest LLC.
789 East Eisenhower Parkway
P.O. Box 1346
Ann Arbor, MI 48106 - 1346
c Copyright by
Vasilis A. Sotiris
2011
Dedication
I dedicate this work to my wife Katya. I owe a lot to my wife, who has
wholeheartedly stood by me, both in good and in bad times, who has helped me,
grounded me, inspired me and believed in me. I would like to thank professor Eric
Slud for his support and guidance throughout my research and for his dedicated
and persistent pursuit to teach and instill mathematical rigor. It’s been an honor to
know and work with professor Slud. I would also like to thank prof. Michael Pecht
standards. Prof. Pecht has taught me the value of knowing the bigger picture,
technical and abstract ideas. I would like to thank Dr. Michael Azarian for always
opening up his time for me and for all his guidance and suggestions. I would like
to thank Dr. Diganta Das for his support, encouragement and always useful insight
and advice. Last but not least, I want to thank all my friends and family for their
support.
ii
Acknowledgments
I would like to acknowledge professor Mei-Ling Ting Lee for her role in intro-
ducing me to first hitting time models and their role in reliability. I would like to
acknowledge Ed Tinsley and Nikhil Vichare at Dell Inc., and Kai Goebel, Abhinav
Saxena and Jose Celaya at NASA, for their support of my research at CALCE.
iii
CONTENTS
1. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Problem Setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Reliability Models using degradation and lifetime data . . . . . . . . 7
1.3 First Hitting Time Degradation Models . . . . . . . . . . . . . . . . . 8
1.4 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.5 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3. ESTIMATION THEORY . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.1 Maximum Likelihood Estimators . . . . . . . . . . . . . . . . . . . . 32
3.2 Methods of Evaluating Estimators . . . . . . . . . . . . . . . . . . . . 33
3.2.1 Finite Sample Measures . . . . . . . . . . . . . . . . . . . . . 33
3.2.2 Asymptotic Evaluations . . . . . . . . . . . . . . . . . . . . . 36
3.3 Asymptotic Properties of Maximum Likelihood Estimators . . . . . . 38
3.4 Computing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.4.1 Observed Fisher Information Matrix . . . . . . . . . . . . . . 39
3.4.2 Multivariate Integration using Gaussian Quadratures . . . . . 39
3.4.3 Nonlinear Optimization . . . . . . . . . . . . . . . . . . . . . 40
v
5.5.1 Contribution to likelihood from failed devices . . . . . . . . . 74
5.5.2 Contribution to likelihood from surviving devices . . . . . . . 77
5.6 General Longitudinal Case - GENL . . . . . . . . . . . . . . . . . . . 79
5.6.1 Contribution to likelihood from failed devices . . . . . . . . . 79
5.6.2 Contribution to likelihood from surviving devices . . . . . . . 81
5.7 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 83
vi
7.5.1 Connection to FHT models . . . . . . . . . . . . . . . . . . . 119
vii
11.5 Multistate Markov Chain . . . . . . . . . . . . . . . . . . . . . . . . . 168
11.6 Estimator for the Transition Probability Matrix . . . . . . . . . . . . 171
11.6.1 Evaluating the Probability Transition Matrix . . . . . . . . . 174
11.6.2 Example of a two-dimensional multistate model . . . . . . . . 175
11.7 Case Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179
11.8 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 186
viii
LIST OF TABLES
x
LIST OF FIGURES
9.1 Algorithm flow diagram showing the processing of the data. . . . . . 132
9.2 Logistic distribution model for posterior class probabilities. . . . . . 142
10.1 Joint posterior class probability vs. observation for Lockheed Martin
test data set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
10.2 Joint posterior class probabilities for CALCEsvm, and the open source
support vector classification software called LibSVM. . . . . . . . . . 147
10.3 Joint positive posterior class probability for simulated data set. . . . 152
10.4 Joint positive posterior class probability for simulated data set in P2. 153
10.5 Joint positive posterior class probability for simulated data set in P3. 154
10.6 Distribution of projected multivariate data on to first two principal
components . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156
10.7 Estimate of unit health as a function of time . . . . . . . . . . . . . . 157
10.8 Variation in estimated health probability across all training units . . 158
10.9 Cross unit variation in the estimated health probability . . . . . . . . 158
10.10GP model fit to the health estimate time series of a test unit . . . . . 159
11.1 Illustration of sample paths from a bivariate stochastic process {X(t), Y (t)}166
11.2 Multistate model representation of a bivariate stochastic process . . . 167
11.3 Shape of probability transition matrix . . . . . . . . . . . . . . . . . 170
xii
11.4 Survival probability estimates from the AJ estimator . . . . . . . . . 178
11.5 Expected time to failure . . . . . . . . . . . . . . . . . . . . . . . . . 179
11.6 Simulation of the degradation process conditioned on an error event
process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
11.7 Expected time to failure distribution for devices with 0,1 and 2 error
events . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183
11.8 Transition probability matrix used in case study . . . . . . . . . . . . 184
11.9 Aalen-Johansen survival probability estimate starting from states 0,
1 and 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185
11.10Expected time to failure starting from state 0,1 and 2 respectively . . 186
xiii
1. INTRODUCTION
The Center for Advanced Life Cycle Engineering (CALCE) at the University of
Maryland has over the past 30 years pioneered methodologies for reliability analysis
of electronic products. With increasing presence of electronics in modern systems
and in every-day products, their reliability is inextricably dependent on that of their
electronics [2]. Reliability methodology for electronics went from simple standard
based assessments, to using physics of failure (PoF) models in the late 1980s, and
more recently prognostics and health management (PHM) models.
The fundamental aim of PoF modeling is to postulate, based on the physics
and mechanics of the failure mechanisms, a set of generic functional relationships
between the mean fatigue life and the operational loads [3]. In their 1990 pa-
per, Dasgupta et al. [3] are among the first to stress the importance of modeling
failure-times in conjunction with PoF models. Pecht et al. 1990 [4], point out the
importance of considering PoF, especially for material properties, for modeling fail-
ures of electronics obtained from laboratory tests. Hu et al. 1991 [5], point out
limitations in conducting accelerated failure tests without knowledge of the failure
mechanisms. They point out the need for ensuring that accelerated failure tests
target and therefore induce the intended failure mechanism.
During use, electronics are exposed to a variety of loading conditions such as
temperature or power excursions, shock and vibration. Interconnects, such as solder
joints, printed circuit board traces, component leads, and connectors are vulnerable
to these loading conditions and are susceptible to failures by mechanisms such as
fatigue, creep, corrosion, and mechanical over-stress [58]. Life cycle loads, either
individually or in various combinations may lead to performance or physical degra-
dation and reduce its service life [7]. Table 1.1 lists life cycle loads experienced by
electronics. The extent and rate of degradation depends on the nature, magnitude,
and duration of exposure to such loads [8].
2
Tab. 1.2: Failure Mechanisms in Electronics
Over Stress Failures Wear-out Failure
Brittle Fracture Wear
Ductile Fracture Corrosion
Yield Dentritic Growth
Buckling Interdiffusion
Large Elastic deformation Fatigue crack propagation
Interfacial De-adhesion Diffusion
Radiation
Fatigue crack initiation
creep
on sample averages. There is also economic and business strategic needs that can be
solved using real-time reliability analysis results. PHM is a method that permits the
assessment of the reliability of a component (or system) under its actual application
conditions [2]. PHM aims to provide advanced warning of failures, enable optimal
maintenance actions, reduce life-cycle costs, and aid in mission critical decisions.
The key element to PHM is its prognostic element, which as we show can play
the role of fault detection, degradation estimation and failure-time prediction. The
measures provided by these predictive outcomes are central to useful implementation
of PHM technology. Pecht 2010 presents a PHM road-map and an assessment of
the state of practice for information and electronic-rich systems [17].
3
Tab. 1.3: Failure modes and mechanisms analysis for the circuit card assembly
Category Site Mode Mechanism Stress
Electrical open thermal Temperature
fatigue cycling
PTH Electrical short Conductive Voltage,
Printed (between PTHs) filament high RH
formation small PTH
circuit spacing
board Electrical short Electro-migration High current
(PCB) (between traces) Corrosion density
Metaliz. and degradation ionic High RH,
traces in resistance contamination electrical bias
open traces
Short between Overheating due
windings Thermal to excessive
and the core fatigue current and
prolonged use at
Compon. Inductors high temperature
Overheating due
Short between Thermal to excessive
windings fatigue current and
prolonged use at
high temperature
Open circuit Thermal Prolonged use
inside the fatigue at high
inductor temperatures
Intermittent Thermal
Inter Solder change in fatigue, Temperature
connect joints electrical creep, high- cycling and
resistance cycle fatigue vibration
the device therefore constitutes a binary process, which is equal to zero while the
device is not failed and equal to one at the failure-time.
This data-structure is deficient in detail and inference models lack predictive
power because of the following reasons:
• A model that only uses lifetime data cannot account for dynamic environments
that products are exposed to in the field. In other words, predictive inference
does not account for real-time environmental or usage conditions.
4
• Often in tests we collect few failures. This can be because we do not have
many devices to test in the first place, typically because they are expensive.
This can also be a result of short test times, constrained by money and time.
Singpurwalla [22] and Sobcyzk [12] among others motivate the use of stochastic
processes in models in order to account for the dynamic environments that fielded
devices experience. In this context, there is therefore a need for data that can
describe the evolution of the stochastic process. In this case PHM prediction models
will become a little more complicated because now they have to accommodate data
measured on the stochastic process. However, these models are presumably better
suited at capturing the changing environments.
Failure tests often result in few failure-time samples. When the failure-time
sample is small, the traditional inference and prediction models suffer because the
available data are not ample enough to fit the model parameters with adequate
precision. There is a need therefore for auxiliary reliability information. Auxiliary
reliability information can come from observing the degradation variable(s) of the
device over time, i.e, the degradation process of the device. The degradation process
consists of a collection of degradation variables indexed by time, and can be thought
of as a stochastic process of accumulating damage. Because failures in electronics
can be strongly linked to known degradation variables (failure modes), and because
the response to stress (even under constant stress) is random (due to variations in
material properties), it is not unreasonable to model degradation as a stochastic
process, and failure as its first hitting time of a threshold.
Due to ”self-healing” of materials in electronics, the damage is not generally
considered non-decreasing over time, but instead, it can ”heal” or recover. Under
5
these physical conditions, a Gaussian process can be appropriate for modeling the
fluctuations in the degradation process. If we can observe this process then we
have access to valuable ”auxiliary” information that can presumably explain the
rate of degradation and therefore improve inference on failure-times. Auxiliary
reliability information can also come from covariates collected on each device, and
when degradation is latent, as we discuss, it can also come from marker variables.
6
So we see that the use of stochastic processes in modeling lifetimes is motivated
from i ) an intuitive representation of degradation, ii ) from a data limitation; namely
the lack of large failure-time samples, and iii ) from the uncertainties generated in
using ALT lifetime data. ALTs do help increase the failure-time sample size, however
they also introduce uncertainty that is difficult to handle in models.
In conclusion, there is a need for PHM reliability models that perform ”well”
under small sample sizes, and we have above pointed out three reasons why we
would have small sample sizes:
• Time/money constraints
• Accelerated test conditions are reduced in order to preserve the failure gen-
erating mechanism. Less stress means fewer samples fail in a fixed period of
time.
It has been noted by Chown, Pullum and Whitmore 1994 [13], that reliance
on lifetime data is becoming less and less practical in engineering, and there exists
a pressing need for reliability models that capture the degradation response of a
device over time. Nair 1988 [15] states: ”... Degradation data are a much richer
source of information than time-to-failure data. The lack of statistical methods for
analyzing them prevents users from exploiting this valuable source of information
[21].
7
Degradation models are based on lifetimes, degradation and covariate mea-
surements that can be collected in the same failure-test. Failure-tests are usually
performed over a fixed time period and some devices may survive. In fact as dis-
cussed in the previous section, it is more common that most devices do survive. The
information on covariates, degradation and lifetimes, collected on both surviving and
failed devices can arguably make PHM technology possible.
Definition 1.1 (Direct failures). Direct failures are defined as the time when an
observable degradation variable first violates a fixed and known failure threshold, so
that the terminal level of degradation is the same for all failed devices.
Definition 1.2 (Indirect failures). Indirect failures are defined as the time when a
latent degradation variable first violates an unknown, possibly random failure thresh-
old, so that the terminal level of degradation varies across devices.
8
like IGBTs, or more basic devices such as capacitors, inductors, resistors, diodes and
transistors. These devices have in common simple structures, for which there exists
an understanding of the physics of failure. Because we understand their PoF, we
know how and why they fail, and therefore can select useful degradation variables
with which to define failure-time.
In failure-tests, indirect failures, are commonly used for electronic systems
or products that are composed of many parts that can interact in complicated and
typically unknown ways. Most modern engineered products depend on sophisticated
electronics, which are housed on densely populated boards or chassis. For electronic
products there are generally no suitable PoF models, and therefore no suitable
degradation variables that can define failure. Failure instead is observed as an
external process, typically by observing the performance of the system, and failure
is defined as the lack of performance to some predefined degree. For example, a
computer freezes, a car stalls, onset of heart attack, etc. In each of these examples
there exists a complicated host system, the computer, the car the human body,
where the true health/degradation is unknown.
We are interested in reliability models that are motivated by both direct and
indirect failures. We are interested in direct failure data because at the component
level we can use PoF models to enhance the predictive power of the data-driven
models. We are interested in indirect failures because they are of greater commercial
and application level importance.
Although direct failures are based on an observable degradation variable, that
variable may not be predictive of failure, i.e., it attains the failure-threshold level
suddenly without any preceding trend. In such cases, the observable degradation
variable is called a surrogate degradation variable, and inference is based on unob-
servable (latent) degradation variables, which we also call true degradation variable.
Some examples of latent degradation variables are: the length of a crack in a solder
9
joint, the surface roughness in a ball bearing fan, etc. These variables are latent
because they cannot be measured during the failure-test.
We are however, interested in the situation (experimental setup) where the la-
tent degradation variable can be measured/determined at failure, potentially through
some intrusive (postmortem or terminal) examination. Access to terminal degra-
dation measurements is important, especially when the test is designed to induce
specific failure mechanisms, and when there are known PoF models to work with. It
is difficult to observe the latent degradation variable during the test without intru-
sive and often destructive procedures, however we can measure it at the failure-time
without interfering. In the example of the solder joint, the crack length can be mea-
sured by cross sectioning and x-ray microscopy analysis, and similarly to determine
the surface roughness on the ball bearing.
Because the latent degradation is only observable at termination, a degrada-
tion model must account for its latency at any other time, motivating what we call
Latent degradation models. We believe that access to true degradation data can help
better estimate the rate and variability of degradation in fielded devices. Using the
true degradation data we gain stronger insight into the effects of the environment
and usage on the rate of degradation. The drawback of using true instead of sur-
rogate degradation data is we only have one such measurement, whereas we have
many observations on the surrogate, on each device. In this case the motivation
for a latent degradation model is only as good as the value/importance of the true
degradation data relative to the surrogate.
For indirect failure data, the need for a latent degradation model is more
obvious. The failure is assumed to occur due to wear-out, but the failure-time is
not determined based on anything we can observe. At the failure-time, however, we
are again interested in the situation where we can measure the terminal degradation
level, on one or more degradation variables.
10
Degradation models should also account for uncertain or unknown failure
thresholds. It is not uncommon in failure-tests to use arbitrary thresholds or thresh-
olds based on antiquated standards. Failure thresholds represent the strength of a
device to sustain stress, and is therefore a heterogeneous quality across devices.
Failure-thresholds may not only vary across devices, but may also vary with time.
In other words, the strength of a device can itself decay/degrade with time.
In failure tests of electronics, reliability data can be observed frequently in time
on each device, and in small failure-time sample situations, it behooves us to use
it. Longitudinal time-indexed observations on device reliability, whether observed
directly on the degradation variable or its surrogate or some higher level marker can
give strong insight on the wear-out and the strength of the device as a function of
time and therefore improve estimation. When the degradation variable is latent, as
we consider here, longitudinal measurements are made on degradation markers that
track the progress of the degradation towards a threshold, giving rise to what we
call bivariate latent degradation models.
11
degradation models. Few have investigated a bivariate stochastic process setting to
model lifetime and degradation data. Predominantly the bivariate structure has
been used to model latent degradation processes, and is found mostly in medical
studies, specifically in immunological and epidemiological studies. There are by far
much fewer literary references to bivariate stochastic processes in the reliability field.
We discuss some of these papers further below.
Sobczyk 1987 [34] and in his 1992 book, presents an exposition of methods
of modeling and analyzing fatigue fracture of engineering materials. He thinks of
fatigue to be random and therefore considers stochastic models for fatigue processes.
He argues that early probabilistic treatment of fatigue was mainly concerned with
the statistics of dispersed data, fitted by various probability distributions, such as
the lognormal and Weibull. He introduces the limitation that such approaches do
not provide any direct relationship to the basic fatigue mechanisms. His model-
ing approach of fatigue consists of three basic steps: (i) choosing an appropriate
stochastic model-process for fatigue accumulation; (ii) determining the probabilistic
properties of the model-process (iii) relating the model process to empirical data
and parameter estimation.
Lu and Meeker 1993 acknowledge small failure samples in engineering failure
tests, specifically in electronics systems. They motivate the need for degradation
models that can be used to define time-to-failure distributions. Based on this idea
they develop statistical methods for using degradation measures to estimate a time-
to-failure distribution for a broad class of degradation models. They use fatigue-
crack-growth data to motivate the model. For each device, they assume that degra-
dation measurements yj are available for pre-specified times tj = t1 , . . . , ts , until
yj crosses the pre-specified critical level D or until a pre-specified censoring time
ts . Sample paths are modeled by a parametric general path model yj = ηj + j ,
j ∼ N (0, σ2 ), and the failure distribution is written in terms of (tj , D, ηj , j ). From
12
here they consider several examples where the model parameters are given specific
parametric forms, such as Weibull, Bernstein, lognormal and multivariate normal.
Kahle 1993 considers the Wiener process as a degradation model of a damage
process. He shows that for independent, not necessary identically distributed ob-
servations of process increments, for observations of lifetime-distributions and for a
mixture of these observations the assumptions of asymptotic normality of Maximum-
Likelihood-Estimations (MLE) are fulfilled. The asymptotic normality of MLEs is
used to find simultaneous confidence regions for parameters of the damage process.
Whitmore 1995 et al. model the degradation process by a Wiener process with
drift. They also model measurement errors as independent normal random outcomes
that are also independent of the degradation process. The true degradation process
is therefore separated from the observed by an error term, that they incorporate
into the model.
In her thesis J. Lu 1995 does an excellent job at motivating the need for
models that use both lifetime and degradation data for estimation and prediction.
She introduces the Wiener process as a degradation model and discusses it for a
mixed data-structure consisting of degradation observations at a set of fixed time
points 0 < t1 < t2 < . . . < tn , and failure-times si . For a failed device i, the mixed
data-structure has the following form: (xi1 , xi2 , . . . , xin , si ), i = 1, . . . , p and for a
surviving device j: (xj1 , xj2 , . . . , xjn ), j = p + 1, p + 2, . . . , p + q. She derives a
likelihood function L(δ, ν), and inference equations for the mixed data-structure.
She also touches on inference when the failure threshold (a) is unknown. Inference
in this case is accomplished by fixing the failure threshold level and estimating
the process parameters for many different such levels. The most suited threshold
level, given the data, can be determined by maximizing the likelihood function
L(δ̂(a), ν̂(a)) over the chosen threshold levels.
In an influential expository paper, Singpurwalla 1995, provides an overview of
13
failure models, based on stochastic processes, that are suitable for describing the
lifetime of items that operate in dynamic [Link] signal a new philosophy
of life-testing experiments wherein one also monitors the environmental factors that
govern tests, and sets a tone for work in the development of models for survival
wherein the physics of failure and the characteristics of the operating environment
play a central role.
Doksum and Normand 1995, present two stochastic models that describe the
relationship between biomarker process values at random time points, event times,
and a vector of covariates. In the first model the biomarker process is a Wiener
process whose drift is a function of the covariate vector. In the second model the
biomarker process is taken to be the difference between a stationary Gaussian process
and a time drift whose drift parameter is a function of the covariates. They present
the methods principally in the context of conducting inference in a population of
HIV infected individuals.
In their 1998 paper, Whitmore et al., present a bivariate Wiener process model
for degradation processes, applied to a terminal data-structure, where degradation is
entirely unobservable. Their paper is one of the earliest sources that use a bivariate
Wiener process to model degradation. This work is the main inspiration for the
model development in part-I of my thesis. Inference is based on observations on a
marker variable that can be used to track the progress of the latent degradation.
They derive joint densities for the likelihood and apply the model to simulated
lifetime and marker data.
Pettit and Young 1999 consider both lifetime and degradation data on failure
and survived devices. They model degradation by a Wiener process, and lifetime as
its first hitting time to a fixed [Link] extend the analysis in J. Lu 1995 by
using a fully Bayesian approach to estimation and prediction.
Lee et al. 2000, extend the bivariate Wiener process considered by Whitmore
14
and co-workers 198, and model the joint process of a marker and a latent health
status. Covariates are related to the model parameters through generalized linear
regression [Link] derive formulas for predicting residual survival time and
discuss model validation on clinical trial data.
Lawless and Crowder 2004 argue that for certain types of degradation processes
a model involving independent non-negative increments is appropriate. They use,
therefore, the Gamma process as a model for degradation processes. They construct
a tractable gamma-process model incorporating a random effect and fit the model
to data on crack growth. Covariates are incorporated via an accelerated life model
by replacing process parameters with a functional of the covariates and a new set
of unknown parameters.
Lee et al. 2006 review first hitting time (FHT) models for survival data, and
introduce threshold regression for survival analysis. They argue that FHT models
can only be valuable in applications if they can include regression structures. Re-
gression structures allow effects of covariates to explain the inherent dispersion of the
data, thereby taking account of variability and sharpening inferences. Threshold-
regression refers to FHT models with regression structures that accommodate co-
variate data. The parameters of the process, threshold and time scales may depend
on the covariates.
Lehman 2008 et al. survey some approaches to model the relationship be-
tween failure time data and covariate data. In particular they consider a class of
degradation-threshold-shock models in which failure is due to the competing causes
of degradation and trauma. They express the failure time in terms of degradation
and covariates, where degradation is modeled by a process with stationary indepen-
dent increments and related to covariates through a random time scale.
Tang and Su 2008, propose to obtain the first hitting times of a degradation
process, modeled by a Wiener process with drift, over certain non-failure thresholds.
15
Based on only these intermediate data, they obtain the uniformly minimum variance
unbiased estimator for the mean lifetime.
Kahle and Lehmann 2010 describe a simple degradation model based on the
Wiener process with drift. They consider the case that each realization of the degra-
dation process, both process increments and failure time are observable, and they
estimate the process parameters. They observe each sample path of the degradation
process to either a failure time or to a censoring time. They develop a likelihood
function based on the conditional distribution of the process under the condition
that the threshold level is not exceeded and the joint distribution of the conditional
process increment and lifetime variable.
Wang and Xu 2010 discuss a class of inverse Gaussian process models for degra-
dation data and associated maximum likelihood inferences. They use an expectation
maximization (EM) algorithm to obtain MLEs of the unknown parameters and the
bootstrap method to assess the variability of the MLEs.
Singpurwalla 2010 provides an interesting perspective on damage accumulation
and marker processes, a perspective and thoughts that are much related to the
ideas developed in this thesis. He talks about damage being an abstract concept,
which is not measurable, but its surrogates can be measured. With this in mind he
highlights a probabilistic architecture based on a bivariate stochastic process with
one component that is non-decreasing and the other that may fluctuate around some
mean. The non-decreasing process leads to the fluctuating observable process. He
argues that the failure threshold is random with an exponential(1) distribution, and
he calls this threshold the hazard potential of an item.
Lee et al. 2010, consider sequential observations on degradation and/or on
covariates prior to failure. They argue there is a need for simple regression methods
to handle longitudinal data, and present the use of the Markov property to do this.
They outline a model that can handle a longitudinal process with an unobservable
16
health status (degradation) as well as time-varying covariates.
A large portion of the literature on degradation models in reliability is found
under accelerated degradation models. The reason as discussed earlier is due to
the need for early failures. Accelerated degradation models also use lifetime and
degradation data collected in ALTs, but differ from ”non-accelerated” degradation
models in that the data structure often includes extra complications that we address
in this work. Work on accelerated degradation models in reliability can safely find
its way back to the middle of the 20th century with early work by Epstein and Sobel
1963, Singpurwalla 1970, 1971 and 1973 [36] [37] [38], Mann et al. 1974 [39] and a
few others, then followed by Bhattacharyya and Fries 1982 [40], Nelson 1990 [41],
Carey and Koenig 1991 [42], Doksum and Hoyland 1992 [116], Meeker and Escobar
1993 [44], Whitmore and Schenkelberg 1997 [45], Lu, Park and Yang 1997 [46],
Meeker, Escobar and Lu 1998 [47], Owen and Padgett 1999 [48] ,Onar and Padgett
2000 [49], Bagdonavicius and Nikulin 2001 and 2004 [50] [51], and more recently,
Padgett and Tomlinson 2004 [52], Park and Padgett 2005 [53], Park and Padgett
2006 [54], Bae, Kuo and Kvam 2007 [55], and Meeker et al. 2009 [56].
Singpurwalla 1970 proposed to investigate the functional relationship between
parameter vector θ for the probability density function of the time-to-failure random
variable and stress vector S. With the above relationship he was interested in mak-
ing inference about the failure behavior of the device at environmental conditions
which cannot be simulated in a test. In this work, he assumes that i) the device
fails due to single failure mode and ii) that the severity of the stress level does not
change the type of life distribution, but that the stress level influences the values of
its parameters. He assumes a linear stress-failure relationship and an exponential
failure time pdf, with the hazard rate parametrized as: λi = BSi , with B unknown
parameter.
Singpurwalla 1973 discusses the problem of inference when both the location
17
and the scale parameter of the time-to-failure distribution are re-parameterized,
the former as a linear function of stress and the latter according to the Arrhenius
re-action rate model. The failure time distribution is again exponential with two
parameters λi and γi , f (t|λi , γi ) = λi exp(t − γi ) and λi = exp(A − B/Vi ), and
γi = α − βVi , where A, B, α and β are unknown parameters. The objective is to
predict the mean time-to-failure at use conditions µu = γu + λ−1
u .
In their book Mann, Schafer and Singpurwalla 1974 discuss accelerated life
testing, models and some results. They are interested in making inference from
accelerated life tests when certain relationships between parameters of a failure
time distribution and the environmental conditions can be reasonably hypothesized.
These relationships or models are derived from an understanding of the physics of
failure (PoF) of the device under discussion. The time-to-failure random variable is
given by f (t|θ) where θ = g(S|a, b, . . .) is known except for a, b, . . . and is valid for
certain ranges of stress S. Their objective is to obtain estimates of a, b, . . . based on
life tests conducted at elevated values/levels of stress, and then use these estimates
to make inference about θ in use environment stress S u .
Bhattacharyya and Fries 1982 focus on the inverse Gaussian IG(θ, λ) as the
failure time distribution. They motivate that the genesis of the IG can be cast
in the context of cumulative fatigue, or depletion of strength. They point out
the relationship to a Wiener process crossing a fixed threshold ω, θ = ωµ−1 and
λ = ω 2 δ −2 , where µ and δ are the Wiener process with drift parameters. In their
accelerated degradation model, the mean of the Wiener process is parameterized as
a linear function of stress as µ = α + βx, ensuring a direct relation between the
cumulative fatigue/wear/degradation and stress levels. For inference they observe
stress and failure time pairs across a range of stress settings.
Carey and Koenig 1991 describe an analysis strategy to extract reliability infor-
mation from measured degradation of devices submitted to elevated stress. Degra-
18
dation is the propagation delay in an integrated logic family device, which increases
with age and temperature. The degradation model on a single device from a given
√
temperature group is: yn − y0 = θ(1 − exp(− λtn )) + n , where n ∼ N (0, σ 2 ),
θ is related to the concentration of impurities in the device and therefore to the
maximum change in propagation delay, yn − y0 is the change observed in the prop-
agation delay between times tn and t0 . The effect of temperature on the maximum
degradation θ is given by log(θ) = A − (B/kT ) + η, where η ∼ N (0, ση2 ), Ti is the
absolute temperature, k the Boltzmans constant and h a random effect representing
unobserved variability.
Doksum and Hoyland 1992 consider step stress accelerated testing, where fail-
ure is modeled in terms of accumulated decay reaching a threshold ω. Accumulated
decay is assumed governed by a Wiener Process W (y). The distribution of W (y)
depends on stress level s(y). The stress level s(y), in turn, is assigned to the device
at each time point y. Time-to-failure is given by Y = IG(y|µ, λ). Their accelerated
degradation model is given by: W (y) = W0 (t + α[y − t]) if y ≥ t and W (y) = W0 if
y < t, where W0 is the Wiener process under the nominal stress level 0. This model
has two stress levels and a decay rate changing from η to αη as y crosses the stress
change point at time t. In the two stress level case, the distribution of the failure
time is IG(τ (y)|µ, λ), τ (y) = y if y ≤ t and equal to t + α(y − t) if y > t. The
model corresponds to making a monotonic transformation of time Y . They think
of Y as the true (calendar) time, and Z = τ (Y ) as the effective (non-accelerated)
time. Lastly in their model they further parameterize α as a function of stress.
Meeker and Escobar 1993 review research and issues in accelerated testing, and
make the point that there are two types of accelerated tests: i) Accelerated Life Tests
(ALT) and ii) Accelerated Degradation Tests (ADTS). In ALTs one observes time-
to-failure information and typically assumes a time-to-failure distribution. In ADTS
one observes at one or more points in time the amount of degradation for a device
19
and typically assumes a model for degradation as a function of time. Traditional
accelerated test statistical models assume a relationship(s) between the constant
stress model parameters.
Whitmore and Schenkelberg 1997 consider a degradation process to be a
Wiener diffusion process with a time scale transformation. The model incorporates
Arrhenius extrapolation for high stress testing. A time transformation accommo-
dates for a time dependent Wiener process drift parameter. The time transforma-
tion depends on the particular degradation mechanism. They encode a relationship
between the parameters and stress level through a functional an Arrhenius model.
Meeker, Escobar and Lu 1998 give a review article on degradation modeling
with accelerated life and degradation data. They assume that the degradation fol-
lows a path defined by: yij = Dij +ij ,i = 1, . . . , n (item) and j = 1, . . . , mi (number
of observations). In this model y is the predicted degradation path, D is the actual
degradation path and ∼ N (0, σ2 ) the error term. They define Dij = h(tij , β i );
β i = (β1i , . . . , βki ), and an Arrhenius model to describe the effect of temperature
on the rate R(temp) = f un(temp) of a simple first order chemical reaction. They
define the acceleration factor AF = R(temp)/R(tempU ). The chemical reaction
model is given as D(ttemp) = D∞ × (1 − exp[−RU × AF (temp) × t]). Solving
for the failure time at temperature T (temp) = T (tempU )/AF (temp). Therefore if
T (rempU ) ∼ W eib(αU , β) then T (temp) ∼ W eib(αU /AF (temp), β).
Owen and Padgett 1999 model the strength of materials with cumulative dam-
age models. They assume that at i) each increment in stress (can be thought of as
time) causes a random amount of non-negative damage D, subject to some distri-
bution function fD and ii) the system has a fixed theoretical strength (threshold)
ψ, but the initial strength W is a random quantity. The cumulative damage after
n + 1 increments in stress is represented by: Xn+1 = Xn + Dn g(Xn ), with g(u) being
the damage model, for example g(u) = 1 gives the additive damage model. They
20
let N be the number of increments of stress applied to the system until failure and
R∞
express the survival probability as: P (N > n) = 0 Fn (w)fW (w)dw. The acceler-
ation variable in their work is the gauge length of the specimen L, knowing that
longer specimens fail faster or equivalently show smaller tensile strength. Therefore,
the distribution of the initial strength variable is parameterized with L.
Onar and Padgett 2001 consider models for the strength of systems based on
cumulative damage arguments. The models are based on a three-parameter inverse
Gaussian distribution, and incorporate system size as a known acceleration variable.
The stress level is denoted as L, the lifetime at stress level L as XL , with cdf FXL (x)
and X ∼ IG(µ, λ). Under a cumulative damage model, a system is placed under
a steadily increasing stress or load until failure occurs. It is assumed that the load
is increased in small discrete increments and that each of these increments causes a
random amount of non-negative damage D > 0, with cdf FD . The systems initial
strength is given by y, and the initial damage by X0 (due to existing flaws or
other damages). The cumulative damage is given as: Xn+1 = Xn + Dn+1 h(Xn ),
and consider h(·) = 1. Again the survival probability is expressed as a function
R∞
of the initial strength W , P (N > n) = 0 Fn (w)fW (w)dw. The initial damage
represents the reduction in strength (or reduction in lifetime) due to inherent flaws
in the system. They assume that the flaw process can be described by a stationary
Gaussian process which given the theoretical strength ψ, yields truncated Gaussian
distribution function for W . They let S represent the total load after N increments
(or calendar time), then S ∼ IG(Λ(θ; L)/ζ, Λ(θ; L)/σ 2 ).
Bagdonavicius and Nikulin (20012) consider degradation models influenced by
covariates under accelerated test conditions. They model degradation by a gamma
process Z(t) = σ 2 γ(t), where γ(t) ∼ Gamma(1, ν(t)) = Gamma(1, m(t)/σ 2 ). They
consider functional forms for the mean degradation m(t) similar to Koenig and
Carey (1991), where m(t) is parameterized in some way m(t) = m(t, g), where
21
g = (g0 , . . . , gm )T are unknown parameters. They assume that the process has
independent increments and therefore the moments are known. A failure caused
by degradation occurs when Z(t) reaches the value z0 , T = inf (t|Z(t) ≥ z0 ).
The stochastic process under the influence of covariates x, is given by: Zx (t) =
Rt
σ 2 γ 0 exp(β T x(s))ds.
In addition to degradation models for latent degradation processes, there is also
a need as we motivate in the thesis, for degradation models on partially observable
degradation processes. In this context, there is no visible literature, and the work
in this thesis we hope can help shed some light and inspire future research.
Bivariate degradation models have predominantly been used to analyze termi-
nal data observations, mostly because degradation is latent. For longitudinal data,
the bivariate model has not been extensively studied. Longitudinal treatment of
markers in the context of bivariate Wiener models can be found in Sy et al. 1997,
Henderson et al. 2000, and Guo and Carlin 2004, among a few selected others,
and predominantly in HIV AIDS studies. No visible literature exists on bivariate
longitudinal degradation models in reliability.
Failure thresholds in most of the work in bivariate degradation models are
considered known and fixed for a given sample of devices. It is however desirable,
and indeed a relevant topic in reliability analysis today, to accommodate random or
uncertain failure thresholds. Literature on random failure thresholds within Wiener
process degradation models, can be found in J. Lu 1995, Singpurwalla, 2010, Wang
and Coit 2007. For bivariate degradation models, however, this is still an open area
of research.
1.5 Contributions
We develop a new degradation model that extends the bivariate Wiener process
model introduced by Whitmore et al. 1998 and allows us to incorporate:
22
1. Terminal degradation observations
23
2. DATA STRUCTURE AND NOTATION
Direct failures are characteristic of electronics, in that although they are de-
fined by an observable degradation variable, that variable is not predictive of failure.
Therefore, through FMMEA and experimentation, more valuable degradation vari-
ables are determined. We consider in the following material, three types of failures
in electronic components:
25
According to Kwon et al. [58], however, DC resistance, often responds too lit-
tle (for example, changes in DC resistance are often obscured by the environmental
noise in a real life situation) or too late (for example, after the crack is large enough
to result in a DC open circuit). Therefore, they argue, DC resistance measurements
would not be expected to provide early indications of interconnect degradation.
Instead, they propose using RF impedance as the degradation variable. They ar-
gue that due to the skin effect, a phenomenon wherein signal propagation at high
frequencies is concentrated near the surface of a conductor, RF impedance exhibits
increased sensitivity to small cracks initiated at the surface of an interconnect. Time
domain reflectrometry (TDR) is a method used to measure RF impedance.
2. A dielectric (or insulator) is a material that resists the flow of electric charge,
and they fail when they collapse (typically due to high voltage) and start to conduct.
The collapse of insulator is sometimes caused by the growth of protruding material
called whiskers, that can grow to create a conductive path between two differently
biased conductors. Tin whiskers are electrically conductive, crystalline structures of
tin that sometimes grow from surfaces where tin (especially electroplated tin) is used
as a final finish. Electronic system failures may occur due to short circuits caused
by tin whiskers that bridge closely-spaced circuit elements maintained at different
electrical potentials. Failure is defined by a DC open circuit, and the degradation
variable used is again DC resistance.
DC resistance, is again not useful for failure prediction since it suffers from the
same effects as in the previous example. Instead, degradation is explained in terms
of the growth of whiskers, for example, the average length of all whiskers larger than
a predefined nominal length (threshold), or the number of whiskers larger than a
predefined threshold, etc. The density and length of whiskers can be measured in a
laboratory setting [59].
26
Fig. 2.2: Whisker growth, seen at initial stage
27
Fig. 2.3: Pre and post aging of an IGBT part. Increased reflectivity picked up by SAM
indicates degradation
The above examples motivate the need for degradation variables better suited
for failure-time prediction. In all three cases failure analysis methods are used to
measure more useful degradation variables, which theoretically are more faithful to
the underlying failure mechanism(s). Some failure analysis methods, like SAM, can
only be performed after failure has been determined, while others, like TDR can be
performed more frequently. In this thesis, we consider only terminal measurements
on degradation and longitudinal observations on markers.
Computers, observed from the system level, which includes both hardware
and software functionality, exhibit indirect failures. Examples of indirect failures in
computers are: sudden shut-down or freeze (blue screen). Hardware variables, like
the motherboard temperature, %CPU throttle, the fan speed, and many more are
28
readily measurable in most personal computers. Software variables are also readily
available, specifically, event indicator variables that flag the occurrence of various
types of programmed errors and warnings. Event processes can be used to model
either the degradation or a marker to degradation over time.
Observed from the system level, Gas-Turbine engines, also exhibit indirect
failures, like when the engine unexpectedly stops producing power. Before engines
stop producing power, they exhibit other intermediate events, that are typically
used as failures. One such event is called compressor surge or stall, which results
when the compressor can no longer compress the incoming air. Typically most
turbine engine failures result from blade degradation due to fatigue and creep of the
materials.
2.2 Data-Structures
The TM data structure is based on Whitmore et al. 1998, and forms the basic
data-structure to which we augment terminal degradation and longitudinal marker
data. We consider one and two marker observations as separate cases. Table 2.1
summarizes the variables used under each data-structure.
29
Tab. 2.1: Summary of variables, their names, description and realizations, under each
data-structure
Terminal Longitudinal
Variable Structure Structure
Name Description TM TMD TMDL
T Event time s∧τ s∧τ s∧τ
=S∧C
S Failure-time s s s
C Censored-time τ τ τ
∆ Failure indicator (0,1) (0,1) (0,1)
=I[S≤C]
Y Marker y(T ) y(T ) y = (y(t0 ), . . . , y(T ))
Z Covariate z(T ) z(T ) z = (z(t0 ), . . . , z(T ))
X Degradation x(s) x(T ) x(T )
A cohort of n independent devices are monitored over a fixed time period [0, τ ],
where the end of testing at τ is considered nonrandom, and known, and typically
chosen as a cost and time constraint. The number of devices that fail by time τ is
P
denoted by q = i=1,...,n I(si < τ ), a binomial random variable that describes the
natural proportion of failed to survived devices in a fixed test period [0, τ ]. The
number of devices that survive by time τ is denoted by p, such that p + q = n. The
failure threshold is again represented by a scalar a, and is assumed known and fixed.
Under the TM data structure, for each device we observe the lifetime T =
min(S, τ ), the marker Y (T ) and the degradation X(T ) at T . For failed devices
we observe T = s, Y (s) and X(s) = a. For surviving devices we observe T = τ ,
and Y (τ ). The TM data-structure is a subset of the data available under the TMD
data-structure. In the TMD data structure we also observe terminal degradation
on surviving devices, that is X(τ ). In both TM and TMD we assume that the
degradation process starts at time t = 0 is equal to zero, X(0) = 0.
30
2.2.2 Longitudinal Structure
31
3. ESTIMATION THEORY
In many cases, there is an obvious candidate for a point estimator, for example,
a sample mean is a natural estimator of a population mean. However, when we leave
a simple case like this, intuition will often lead us astray. Therefore, it is useful to
have techniques of arriving at reasonable candidates for consideration.
The method of maximum likelihood is by far the most popular technique for
deriving estimators. The following are some of the advantages of this estimation
technique. Maximum likelihood provides a consistent approach to parameter esti-
mation problems. This means that maximum likelihood estimates can be developed
for a large variety of estimation situations. For example, they can be applied in
reliability analysis to censored data under various censoring models. Also, it is well
known that maximum likelihood methods have desirable mathematical and opti-
mality properties. We will discuss these properties in the next section.
The disadvantages of maximum likelihood estimation include the following.
The likelihood equations need to be specifically worked out for a given distribution
and estimation problem. The mathematics is often non-trivial, particularly if confi-
dence intervals for the parameters are desired. The numerical estimation is usually
non-trivial. Maximum likelihood estimates can be heavily biased for small samples
and can be sensitive to the choice of starting values. Maximum likelihood estimation
requires the adoption of strong assumptions about the joint density of the data, if
the assumptions fails, MLEs may be inconsistent.
Consider the likelihood function L(θ|X) = L(θ1 , . . . , θk |X1 , . . . , Xn ).
For ease of exposition, in this section we assume that the parameter vector
consists of a single parameter θ. All of the results presented here extend to the case
of a multi-parameter distribution.
Definition 3.3 (Mean Squared Error). The mean squared error (MSE) of an
estimator W of a scalar parameter θ is the function of θ defined by: Eθ (W − θ)2 .
33
Other distance measures between the parameter and its estimator can be used
as a measure of performance of an estimator. In general, any increasing function
of the absolute difference |W − θ| can be considered, but the MSE has at least two
advantages over other distance measures. First, it is quite tractable analytically
and, second, it has the interpretation:
Thus, the MSE incorporates two components, one measuring the variability
of an estimator and the other its bias (accuracy). To find an estimator with good
MSE properties, we need to find estimators that control both variance and bias.
Z
d ∂
Eθ W (X) = [W (x)f (x|θ)]dx (3.1)
dθ ∂θ
34
d 2
Eθ W (X)
V arθ (W (X)) ≥ dθ (3.2)
∂ 2
Eθ logf (X|θ)
∂θ
This theorem specifies the lower bound on V arθ W , thus an estimator W which
attains this variance is UMVUE. Such an estimator is also referred to as finite-
sample efficient. The quantity Eθ ((∂/∂θlogf (X|θ))2 ) is called the Information
number or Fisher information of the sample. As the information number gets
bigger and we have more information about θ, we have a smaller bound on the
variance of the best unbiased estimator.
∂ ∂
I(θ) i,j
=E logf (X|θ) logf (X|θ) θ (3.3)
∂θi ∂θj
The Cramer-Rao lower bound for θi is then the ith diagonal element of the
inverse of the Information Matrix.
35
Lemma 3.1. If f (x|θ) satisfies
Z
d ∂ ∂ ∂
Eθ logf (X|θ) = logf (x|θ) f (x|θ) dx
dθ ∂θ ∂θ ∂θ
then
∂2
2
Eθ ∂/∂θlogf (X|θ) = −Eθ logf (X|θ)
∂θ2
Informally, this says that as the sample size becomes infinite, the estimator
will be arbitrarily close to the parameter with high probability. Equivalently, we
can say that the probability that a consistent sequence of estimators misses the true
parameter is arbitrarily small (or converges to 0).
36
√ D
ing in proportion to 1/ n as the sample size n grows. Using −
→ to denote conver-
gence in distribution, Wn is an asymptotically normal sequence of estimators if
√ D
n(Wn − θ) −
→ N (0, V )
In practice, we do not deal with infinite samples and therefore it makes sense
to talk instead of ”large-sample efficiency”. Informally, an estimator Wn of θ is
large-sample efficient if, for a large enough sample size n, its empirical distribution
is centered around the true value of θ and it’s empirical variance is equal to V /n plus
a remainder of smaller order than 1/n, where V is the approximated Cramer-Rao
lower bound. To approximate the Cramer Rao lower bound for a given sample size
n, we approximate the expected information number with the observed information
number, −∂ 2 /∂θ2 logL(θ|X)|θ=θ̂ .
Large-sample efficiency of a sequence of estimators Wn (X1 , . . . , Xn ) can be
checked via simulation by sampling k independent vectors (X1n , . . . , Xkn ) of size n,
computing Wn = (Wn1 , . . . , Wnk ) and studying the sampling distribution of Wn as n
grows. For an asymptotically efficient sequence, we expect the sampling distribu-
tion to look more and more normal, centered around the true parameter, with the
empirical estimate of the variance asymptotic to V /n where V is the approximated
(via observed information number) Cramer Rao lower bound for large n.
37
Definition 3.8 (Asymptotic Relative Efficiency). If two estimators Wn and Vn
√ 2
√
satisfy n[Wn − g(θ)] → N [0, σW ], and n[Vn − g(θ)] → N [0, σV2 ] in distribution,
the asymptotic relative efficiency (ARE) of Vn with respect to Wn is:
2
σW
ARE(Vn , Wn ) = (3.4)
σV2
regularity conditions [61] on f (x|θ) and, hence, L(θ|x), for every > 0 and every
θ ∈ Θ,
limn→∞ Pθ (|θ̂ − θ| ≥ ) = 0
where V (θ) is the Cramer-Rao Lower Bound. That is, θ̂ is a consistent and asymp-
totically efficient estimator of θ.
38
Note that asymptotic efficiency is defined only when the estimator is asymp-
totically normal and, asymptotic normality implies consistency.
3.4 Computing
39
We calculate nodes x1 , . . . , xn and coefficients w1 , . . . , wn such that
Z b N
X
f (x)dx ≈ wi f (xi ) (3.7)
a i=1
Given the nodes xi , i = 1, . . . , N one can find the weights wi by solving a set
of linear equations for the weights such that the quadrature (3.7) gives the correct
answer for the integral of the first N − 1 orthogonal polynomials. There are more
efficient ways of computing weighting coefficients, such as through the eigenvalue
decomposition of the symmetric, tridiagonal Jacobi matrix. For further details we
point to the following references [62].
We implement Gaussian quadratures by first selecting N , and computing xi
and wi , i = 1, . . . , N , and storing the vectors into memory. When integration is
required we apply (3.7) using the stored nodes and corresponding nodes.
40
4. WIENER PROCESS AS A DEGRADATION MODEL
The joint densities in the Whitmore and our extended model are of mixed
type. They require the joint relationship between continuous and discrete random
variables. Degradation and marker variables are considered continuous type ran-
dom variables, while the failure-indicator variable is considered a discrete random
variable. In this section we define some of the machinery for specifying mixed type
joint densities.
Let U be a continuous random variable with density f (u) and let V be a
discrete random variable taking values vi , i = 1, 2, . . . , n, with probabilities P (vi ).
To characterize the relationship between U and V we specify the conditional density
f (u|vi ) or the conditional probability P (vi |u). For any a ≤ b the conditional density
f (u|vi ) must satisfy:
Z b
f (u|vi )dx = P (a ≤ U ≤ b|V = vi )
a
Assuming that all continuous-variable densities and joint densities are contin-
uous functions of their arguments u, the conditional probability P (vi |u) is defined
as the limit of P (vi |u ≤ U ≤ u + ) as goes to zero.
If f (u) and f (u|vi ) are all continuous function of u, then by the mean value
theorem:
f (u|vi )P (vi )
P (vi |u) = (4.1)
f (u)
42
(1, 2, . . . , n) and V = {vi |i ∈ I}.
RbP P
a i∈I f (u|vi )P (vi )dx = i∈I P (a ≤ U ≤ b|V = vi )P (V = vi )
P
= i∈I P (a ≤ U ≤ b, V = vi ) = P (a ≤ U ≤ b, V ∈ V)
The product f (u|vi )P (vi ) plays the role of the bivariate function ψ(u, vi ).
Hence by equation (4.1) so does P (vi |u)f (u). Material on mixed-type densities
is taken from [71].
The basic degradation model used in this work is given by the Wiener process
with drift:
X(t) = x0 + νX t + σX WX (t) (4.2)
2
A Wiener process X(t) with drift νX and variance σX has stationary and
independent Gaussian increments with probability density given by:
1 x − x 0 − νX t
fX(t) (x) = √ φ √ (4.3)
σX t σX t
It is assumed that each device experiences its own degradation process which
is independent of other devices. Devices from the same ”family”, having the same
design, are assumed to have the same drift and variance parameters. Covariates can
be used to model the heterogeneous drift in the degradation process as a function
of some other parameters α, β, . . .. In electronics, accelerated lifetimes are most
common, induced by higher than normal experimental stress factors z, such as
temperature, humidity and pressure. It becomes important to model the influence
of such stressors on the rate of degradation. We discuss this further in chapter 7.
43
4.3 Definition of a Wiener Process
[63]
• Random variables Wt0 , Wt1 − Wt0 , . . . , Wtk − Wtk−1 are independent for every
k ≥ 1 and 0 = t0 ≤ t1 ≤ . . . ≤ tk .
[63]
44
Remark 4.1 (Distribution). Since Wt − Ws ∼ N (0, |t − s|), it follows that ξ :=
(Wt − Ws )|t − s|−1/2 ∼ N (0, 1) for t 6= s and the distribution of Wt − Ws has the
density:
x2
1
f (x) = p exp − (4.4)
2π|t − s| 2|t − s|
k
! k
\ Y
Pr Ai = P r(Ai )
i=1 i=1
for all Ai ∈ Fi .
According to lemma 4.1, there always exists a filtration with respect to which
Wt is a Wiener process.
45
Definition 4.6. Let Wt be a {Ft }-adapted Wiener process; assume that for any t,
h ≥ 0 the random vector Wt+h − Wt and σ-algebra Ft are independent. Then we will
say that Wt is a Wiener process with respect to {Ft }, or that (Wt , Ft ) is a Wiener
process.
Theorem 4.2 (The Markov Property). Let (Wt , Ft ) be a Wiener process. Fix
t, h1 , . . . , hn ≥ 0. Then the vector (Wt+h1 − Wt , . . . , Wt+hn − Wt ) and the σ-algebra
Ft are independent. Furthermore, Wt+s − Wt , s ≥ 0, is a Wiener process [67].
Theorem 4.2 says that, for every fixed time t ≥ 0, the process Wt+s − Wt starts
fresh as a Wiener process ”forgetting” everything that happened before time t. For
a Wiener process, this property has a natural extension when t is replaced with
a random time s, provided that s does not depend on the future in any way. To
describe exactly what we mean by this, we continue with the following definition:
Theorem 4.3 (Strong Markov Property). Let (Wt , Ft ) be a Wiener process and
s an Ft stopping time. Assume that P (s < ∞) = 1. Let
W
F≤s = σ{{ω : Wu∧s ∈ B}, u ≥ 0, B ∈ B}
W
F≥s = σ{{ω : Ws+u − Ws ∈ B}, u ≥ 0, B ∈ B}
W W
Then the σ-algebras F≤s and F≥s are independent in the sense that for every A ∈
W W
F≤s and B ∈ F≥s we have P (AB) = P (A)P (B). Furthermore, Ws+t − Ws is a
Wiener process [67].
The strong Markov property gives justification to one of the most important
properties of Wiener processes, a property that forms the basis for the likelihood
46
equations to follow; the Reflection Principle. The Reflection principle helps simplify
otherwise complicated probability expressions related to the Wiener process.
Proposition 4.1 (The Reflection Principle for a Wiener process with zero
drift). Let W (t) be a Wiener process with νX = 0, a > 0, and sa = inf {t : W (t) ≥
a}, then:
fW (t),I(S<t) (w, 1) = fW (t) (2a − w)
The argument for proposition 4.1 is made by noticing that if sa < t, then W (t)
is conditionally just as likely to be above or below level a by the same distance.
Proposition 4.2 (The Reflection Principle for Wiener process with drift).
Let X(t) be a Wiener process with X(0) = 0, νX 6= 0, a > 0, and sa = inf {t :
X(t) ≥ a}, then:
2νX (x − a)
fX(t),I(S<t) (x, 1) = exp 2
fX(t) (2a − x) (4.5)
σX
The joint density in equation (4.5) is a building block for constructing like-
lihood equations under the longitudinal data-structures. We call this density the
complimentary Wiener term, and we derive it in section 4.7.3.
Like the Weibull and logNormal distributions the Inverse Gaussian distribu-
tion is used to model lifetime data. Chhikara and Folks 1977 [68] study the use of
the inverse Gaussian distribution for a lifetime model and suggest its application for
studying reliability aspects when there is a high occurrence of early failures. Unlike
the Weibull and log Normal, the inverse Gaussian distribution is also physically
justified as the first hitting time of a Wiener process to a threshold, which implies a
natural applicability in studying degradation processes which lead to failure events.
47
This means that the lifetime distribution described by the inverse Gaussian density
function depends on the drift and the variance of the degradation process. This re-
lationship facilitates the development of degradation models that use mixed lifetime
and degradation data and provide a natural framework of incorporating covariates.
Definition 4.8 (First Hitting Time). Define Sa , the first time at which the
process X(u), starting from x0 = 0, with drift νX ≥ 0 reaches level a, by: Sa =
inf {u ≥ 0 : X(u) ≥ a}, a > x0 . If a 6= 0, then the distribution of Sa has the density
f (s) = 0, for s ≤ 0, and , ∀ s > 0.
(a − νX s)2
2 −1/2
fSa (s) = (2πσX ) |a|s−3/2 exp − 2
(4.6)
2σX s
The density in equation (4.6) is called the Wald density or the inverse-Gaussian
density. The inverse-Gaussian density can also be expressed in terms of its scale and
shape parameters µ and λ.
1/2
λ(s − µ)2
λ
fSa (s; µ, λ) = exp − (4.7)
2πs3 2µ2 s
48
by measuring variables which characterize the degradation (aging) of each observed
device over time.
In practice, degradation variables in electronics are typically latent (unobserv-
able). For this reason, lifetime predictions are typically based on information from
observed marker variables that track degradation. Markers, unlike covariates are
random variables related to degradation through a parametric model with unknown
parameters. For example, in printed circuit boards, the degradation variable for
a degrading solder joint can be the length of a crack. It is impractical and often
impossible to measure the length of a crack, and we depend therefore on markers
to the crack-length, such as resistance and capacitance across the joint. From a
cost perspective, degradation information is seen as being expensive, and Marker
information as being cheap.
Our basic Marker model is given by a Wiener process with drift. The mo-
tivation for a Wiener marker process is mathematical convenience in expressing a
bivariate marker-degradation relationship. The critical component of the marker-
degradation model is the correlation coefficient, made available through the joint
Gaussian relationship. Our basic marker model, like the degradation model, is
given by a Wiener process with drift:
1 y − y0 − νX t
fY (t) (y) = √ φ √ (4.9)
σY t σY t
Following [24] and [26], the basic analytical framework for a degradation-
marker process is an independent-increments bivariate process, in which all paired
increments X(t) − X(s) and Y (t) − Y (s) have correlation coefficient ρ, which does
49
not vary over time.
The vector {X(t1 ), Y (t2 )} has a bivariate normal distribution with mean vector
(νX t1 , νY t2 )0 and variance covariance matrix ΣXY .
2
t1 σX t1 ∧ t2 ρσX σY
ΣXY = (4.10)
t1 ∧ t2 ρσY σX t2 σY2
We assume that ΣXY is positive definite, νX ≥ 0 and |ρ| > 0 and close to
1. Weak correlation between the marker and the degradation variables will reduce
the predictive efficiency of the model. With c = ρσY /σX , the conditional density of
Y (t) given X(t) is:
!
1 y − y0 − νY t − c(x − νX t)
fY (t)|X(t) (y|x) = p φ p (4.11)
σY (1 − ρ2 )t σY (1 − ρ2 )t
X µX ΣXX ΣXY
∼ N ,
Y µY ΣY X ΣY Y
where
µX ∈ Rm , µY ∈ Rn , and the matrix blocks ΣXX , ΣXY , ΣY X and ΣY Y . are
respectively of size m × m, m × n, n × m, n × n.
50
where, for x0 = 0, y0 = 0 we have:
µY |X = νY r + ΣXY Σ−1
XX (x − νX t)
ΣY |X = ΣY Y − ΣXY Σ−1
XX ΣY X
Within the bivariate-Gaussian case (X(tj ), Y (tj )), Q(tj ) ∼ N (µQ (tj ), ΣQ (tj ))
where µQ (tj ) = tj (νY − cνX ), and ΣQ (tj ) = σY2 tj (1 − ρ2 )
51
Proof. Let Σj,i = Cov(Qj , X i ). With c = ρσY /σX , Σj,i = 0. Therefore, Qj ⊥ X i
due to the fact that a zero covariance implies independence for jointly Gaussian
vectors.
Corollary 4.2.1. The entire process Q(·) is independent of the entire process X(·)
The corollary 4.2.1 of lemma 4.2 is stated in a more general way so that the
finite number of evaluation points for Q are not necessarily the same as the finite
number of evaluation points for X.
Lemma 4.3. Process Q(·) ⊥ (X(·), I(S ≥ t)) where I(S ≥ t) = I(X(u) ≤ a, u ≤ t)
Proof. This follows from corollary 4.2.1 and measurability of Sa with respect to X i
Lemma 4.5. Let S̃a be a shifted first-hitting time defined as S̃a (u) ≡ Sa (u − tj ),
then, fSa |X(tj ),I(S≥tj ) (u|x, 1) = fS̃a−x (u), u ≥ tj , where
!
a−x (a − x − νX (u − tj ))2
fS̃a−x (u) = p exp − 2
(4.13)
2
2πσX (u − tj )3 2σX (u − tj )
Proof.
Sa |(X(tj ) = x, I(S ≥ tj ) = 1) =
= inf {u ≥ tj : X(u) ≥ a|X(tj ) = x, S ≥ tj }
= inf {u ≥ tj : X(u − tj ) ≥ a − x|X(tj ) = x, S ≥ tj } (1)
= inf {v ≥ 0 : X(tj + v) − X(tj ) ≥ a − x|X(tj ) = x, S ≥ tj }
52
X
By theorem 4.3, because X(tj + v) − X(tj ) is F≥tj
measurable, it is independent
X
of (X(tj ), I(S ≥ tj )) which is F≤tj
measurable. Therefore we can simplify (1) to
inf {u ≥ tj : X(u − tj ) ≥ a − x} = inf {u − tj ≥ 0 : X(u − tj ) ≥ a − x} ∼ S̃a−x .
Generalizing Y away from being Wiener, is also of interest, especially when the
relationship between X(·) and Y (·) is non-Gaussian. Such a generalization is rea-
sonable as long as Y (·)−cX(·) is selected to be any independent-increments process,
independent of X(·), with X(·) Wiener with drift, but Y (·) − cX(·) distributionally
unspecified. We propose this extension as part of future research.
53
function is given by:
q
Y
Lθ = fYSa ,XSa ,T,I(Sa <τ ) (Yi (Si ), a, Si , 1)×
i=1
Y p (4.14)
fYτ ,Xτ ,T,I(S<τ ) (Yj (τ ), Xj (τ ), τ, 0)
j=1
For a given realization of the random variables above, the likelihood be-
comes a function of parameter vector θ. For failed devices we observe (Yi (Si ) =
yi , Si = si ; Xi (Si ) = a), i = 1, . . . , q and for surviving devices we observe (Yj (τ ) =
yj , Xj (τ ) = xj ), j = 1, . . . , p.
Tab. 4.1: Conditional-density terms under TMD for ith device, s < τ
Failed devices: fYS ,XS ,T,I(S<τ ) (y, a, s, 1) =
=fYS |XS ,S (y|a, s)fS (s)
=fQs (y − ca)fS (s)
Surviving devices: fYτ ,Xτ ,T,I(S<τ ) (y, x, τ, 0) =
=fQτ (y − cx)fXτ ,I(S<τ ) (x, 0)
where fXτ ,I(S<τ ) (x, 0) = fXτ (x) − fXτ ,I(S<τ ) (x, 1)
and fXτ ,I(S<τ ) (x, 1) = (4.5) by proposition 4.2
In the following sections we derive analytical equations for the joint density
terms required for the likelihood function. Table 4.1 summarizes the required terms
for failed and surviving devices. For short we use the notation S in place of Sa .
54
with the probability density function of S given by equation (4.6), and the proba-
bility density of Qs given by:
!
1 u − y0 − c(a − x0 ) − (νY − cνX )s)
fQs (u) = p φ p (4.16)
σY (1 − ρ2 )s σY (1 − ρ2 )s
Proof. For failed devices, the terminal time random variable T and lifetime random
variable Sa are equal, with values denoted si < τ . Therefore, for s < τ ,
fYS ,XS ,T,I(S<τ ) (y, a, s, 1) = fYS |XS ,S,I(S<τ ) (y|a, s, 1)fXS ,S,I(S<τ ) (a, s, 1) (4.17)
For S = s < τ , I(S < τ ) = 1 is a degenerate random variable, the first term above
can therefore be written as:
The second term in equation (4.17) can be simplified through the observation that
on the event (S < τ ), (XS , I(S < τ )) is a degenerate random-variable pair equal by
definition to (a, 1), and therefore the density fXS ,S,I(S<τ ) (a, s, 1) can be replaced by
fS (s)
fXS ,S,I(S<τ ) (a, s, 1) = fS (s)
55
Finally, the contribution to the likelihood from a failed device is given by:
The contribution to the likelihood from data on surviving devices is given by:
Proof. For surviving devices, the non-degenerate observables are (yj , xj ). Therefore,
variable T = τ can be dropped from the density.
fYτ ,Xτ ,I(S<τ ) (y, x, 0) = fYτ −cXτ +cXτ |Xτ ,I(S<τ ) (y|x, 0)fXτ ,I(S<τ ) (x, 0) =
by lemma 4.3
= fQτ |Xτ ,I(S<τ ) (y − cx|x, 0) = fQτ (y − cx)
therefore,
fYτ ,Xτ ,I(S<τ ) (y, x, 0) = fQτ (y − cx)fXτ ,I(S<τ ) (x, 0) (4.19)
The first density factor on the right-hand side of equation (4.19) is normal,
with s replaced by τ . The second, is the probability density function at x (necessarily
56
< a) of the terminal value X(τ ) for a device surviving at time τ . Any degradation
sample path starting at X(0) = x0 < a and terminating at X(τ ) = x < a either
does or does not cross the failure threshold at a > 0 in the interval (0, τ ). By the
law of total probability:
Equation (4.20) says that the probability of reaching a terminal value x for a
device surviving at time τ is equal to the probability of a Wiener process with drift
reaching a value of x at time τ minus the probability of reaching x and crossing the
threshold at some time earlier.
Using equation (4.5) for the complimentary Wiener term, fXτ ,I(S<τ ) (x, 1), we
get:
2νX (x − a)
fXτ ,I(S<τ ) (x, 0) = fXτ (x) 1 − exp 2
(4.21)
σX
Plugging in equation (4.21) into equation (4.18) we get the final expression for
the likelihood contribution from data on surviving devices:
2νX (x − a)
fYτ ,Xτ ,I(S<τ ) (y, x, 0) = fQτ (y − cx)fXτ (x) 1 − exp 2
(4.22)
σX
57
by:
Z τ
fXτ ,I(S<τ ) (x, 1) = fXτ −Xs (x − a)fS (s)ds (4.23)
0
(x − a − νX (τ − s))2
1
fXτ −Xs (x − a) = p exp − 2
2
2πσX (τ − s) 2σX (τ − s)
(a − νX s)2
a
fS (s) = p exp − 2
2πσX 2 3
s 2σX s
Proof. By definition, fXτ ,I(S<τ ) (x, 1)dx is the probability the degradation level be-
longs to a small interval (x, x+dx) at time τ , and that it crossed the failure threshold
some time earlier.
τ
P (Xτ ∈ (x, x + dx), S = s)
Z
= ds (4.24)
0 dx
X
By theorem 4.3, because Xτ − Xs is F≥s measurable, it is independent of
X
(X(s), S) which is F≤s measurable. Therefore, we can simplify as follows:
Plugging in equation (4.26) into (4.24) we get the final analytical expression
58
for the complimentary Wiener density term:
Z τ
fXτ ,I(S<τ ) (x, 1) = fXτ −Xs (x − a)fS (s)ds (4.27)
0
2νX (x − a)
fXτ ,I(S<τ ) (x, 1) = exp 2
fXτ (2a − x) (4.28)
σX
Proof.
Z τ
fXτ ,I(S≤τ ) (x, 1) = fXτ |S (x|s; νX )fS (s; a)ds (4.29)
0
(x − a − νX (τ − s))2
= fXτ −Xs −νX (τ −s) (x − a − νX (τ − s); 0) = exp − 2
2σX (τ − s)
We take the ratio of the last expression at νX over the same expression with
νX replaced by 0.
(x − a − νX (τ − s))2
exp − 2
2
2σX (τ − s) 2νX (x − a) − νX (τ − s)
C1 (x, s) = = exp
(x − a)2 2
2σX
exp − 2
2σX (τ − s)
Therefore,
2
2νX (x − a) − νX (τ − s)
fXτ |S (x|s; νX ) = exp 2
fXτ |S (x|s; 0)
2σX
Above, by theorem 4.3, fXτ |S (x|s; 0) = fXτ −Xs (x − a; 0), which by proposition
59
4.1 is equal to fXτ |S (2a − x|s; 0). Then,
2
2νX (x − a) − νX (τ − s)
fXτ |S (x|s; νX ) = exp 2
fXτ |S (2a − x|s; 0) (4.30)
2σX
where
(a − x)2
exp − 2
fXτ |S (2a − x|s; 0) 2σX (τ − s)
C2 (x, s) = =
(a − x − νX (τ − s))2
fXτ |S (x|s; νX )
exp − 2
2σX (τ − s)
Plugging in C2 into (4.30) and the resulting (4.30) into (4.29) by direct calculations
we get:
2νX (x − a)
fXτ |S (x|s; νX ) = exp 2
fXτ |S (2a − x|s; νX ) (4.32)
σX
Z τ
2νX (x − a)
fXτ ,I(S<τ ) (x, 1) = exp 2
fXτ |S (2a − x|s; νX )fS (s)ds (4.33)
σX 0
The integral in equation (4.33) is equal to: fXτ ,I(S<τ ) (2a − x, 1). Since 2a − x > a
we get finally the expression in equation (4.5):
2νX (x − a)
fXτ ,I(S<τ ) (x, 1) = exp 2
fX(τ ) (2a − x) (4.34)
σX
60
a Wiener process path which is reflected symmetrically around the level a after the
instant S of hitting a.
The likelihood function in equation (4.14) is computed with factors given by
equation (4.15) for devices failing before τ , and by equation (4.18) for devices sur-
viving past τ .
We aim to predict the degradation level and failure-time for a device surviving
at time t, and whose marker and covariate vector are known at time t. The condi-
tional density of the degradation variable at time t and of the failure-time density
at future time s ≥ t are given by:
and
g(s|y; θ) = fS|Yt ,I(S<t) (s|y, 0) (4.36)
Related, are more complex expressions for density functions of X given a vec-
tor observation on Y . Some of this material is discussed in chapter 5. For now
we derive analytical expressions for the above density functions given one marker
observation. Predictive inferences are based on MLEs of the process parameters,
and should therefore consider sampling error and predictive uncertainty. Neverthe-
less, useful insights are obtained by considering predictive inference when process
parameters are assumed known. By varying the response variable across its domain
we can compute the density at each discrete sample point. For example, we evaluate
function g(s|y; θ), for a given θ on the range 0 ≤ s ≤ ∞ to get an understanding of
the distribution of g(·).
61
Lemma 4.8 (Predicted degradation density). The probability density function
of the degradation random variable conditioned on the marker and survival at time
t is given by:
2νX (x − a)
fQt (y − cx)fXt (x) 1 − exp 2
σX t
h(x|y; θ) = (4.37)
Ra 2νX (u − a)
f (y − cu)fXt (u) 1 − exp
−∞ Qt 2
du
σX t
fXt ,Yt ,I(S≥t) (x, y, 1) = fYt |Xt ,I(S≥t) (y|x, 1)fXt ,I(S≥t) (x, 1) (4.39)
The first factor is given by equation (4.11) due to lemma 4.4, and the second
factor by equation (4.21). The denominator in equation (4.38) is the joint density
of the marker and survival at time t, and is computed by integrating out X(t).
Z a
fYt ,I(S≥t) (y, 1) = fY (t)|X(t),I(S≥t) (y|x, 1)fXt ,I(S≥t) (x, 1)dx (4.40)
−∞
then we have
2a(x − a)
fQt (y − cx)fXt (x) 1 − exp 2
σX t
h(x|y; θ) = (4.41)
Ra 2a(u − a)
f (y − cu)fXt (u) 1 − exp
−∞ Qt 2
du
σX t
62
of the future failure-time random variable, conditioned on the marker and survival
at time t is given by:
Ra 2a(u − a)
−∞ S
f (s − t; a − u)fQt (y − cu)fXt (u) 1 − exp 2
du
σX t
g(s|y; θ) =
Ra 2a(u − a)
f (y − cu)fXt (u) 1 − exp
−∞ Qt 2
du
σX t
! (4.42)
2
a−x (a − x − νX (s − t))
fS (s − t; a − x) = p exp − 2
2πσX2
(s − t)3 2σX (s − t)
Z a
fS,Yt ,Xt ,I(S≥t) (s, y, x, 1)dx =
−∞
Z a
= fS|Xt ,Yt ,I(S≥t) (s|x, y, 1)fXt ,Yt ,I(S≥t) (x, y, 1)dx (4.44)
−∞
63
5. LONGITUDINAL MARKER MODEL
1. The presence of repeated measurements for each device implies that the ob-
servations from the same device are autocorrelated or serially correlated. This
requires us to develop statistical methodology that takes the serial correlation
into account.
3. Most longitudinal data from practical studies contain missing data. Dealing
with missing data, when the missing data mechanism is informative, is gen-
erally nontrivial. To make proper statistical inference, one has to rely on the
information that is supposed to be, but actually not, observed.
65
grating out the latent variable. With the availability of the joint distribution, the
full maximum likelihood estimation and inference can be developed. Albert (1999)
points out that this modeling approach is particularly useful for analyzing longitu-
dinal data in which there is a sizable number of missing observations either due to
missed visits, loss to follow-up, or death. Conditional modeling via a latent variable
can pose some challenges:
1. If the dimension of the latent variable is high, numerical evaluation of the joint
distribution in the likelihood function can be intricate, which typically leads
to computational problems.
In this data setting, time plays a more important role in the likelihood func-
tion and contributes therefore more strongly to parameter estimation. With highly
reliable manufactured products, like modern electronic components, for example,
very few failures are observed in tests. On the other hand, rich longitudinal covari-
ate histories are typically collected and stored during the test. Longitudinal marker
models, are therefore becoming increasingly more relevant and in demand in the
area of reliability.
The intuition and hypothesis is that additional information on degradation
markers will improve the accuracy of parametric and predictive inference as com-
pared to inference in the terminal data-structure case, discussed in chapter 4. We
derive the joint densities for failed and surviving devices under TMDL, just as in
the TMD data-structure. In this case, however, the joint density must account for
the dependency between marker observations in time. After all, the value of the
marker at a certain time point is certainly dependent on the history of the marker
process up until that time.
66
5.2 Scheduling Longitudinal Measurements
Typically marker variables are monitored using appropriate sensors, that can
collect data at very high frequencies. When marker observations are cost or time
prohibitive, failure-tests must be designed accordingly. One of the natural questions
that arises is: ”How many marker observations are needed for this test?”. This
question assumes there exists a minimum number of necessary marker observations,
after which any additional marker observations do not improve parametric inference.
This is an important question in many fields. From a medical stand point,
in designing clinical trials, for example, practitioners may ask patients to make
scheduled return visits to the hospital in order to collect marker data over a certain
period of time. Complications arise when patients don’t return on schedule, or miss
their scheduled appointments. In addition, data on markers can require complicated,
intrusive and expensive procedures. It is therefore of interest to intelligently plan
for the minimum number of scheduled visits. In reliability tests the idea is very
similar. Data on degradation markers can also require complicated, intrusive and
expensive measurement procedures, and examples abound in various applications.
In electronics, however, although marker measurements are typically easy to collect,
there is still interest in reducing the size of the collected data in order to expedite
computations.
Another important question is when to schedule marker measurements. For
a given device, is it better to collect data on the marker(s) with a fixed frequency,
or given a fixed number of possible observations, is it better to observe the marker
more frequently later or earlier in the device’s life?
67
5.3 Parametric Inference
We show that expressions for the joint density between marker observations
have the same form for both failed and survived items. What differs is only the
last term, which accounts for the lifetime observation T = min(S, τ ). Table 5.1
summarizes the key density terms used in the likelihood function for the TMDL1
data-structure.
68
5.4 One Intermediate Marker Observation - TMDL1
q
Y
Lθ = fYS ,Y1 ,XS ,T,I(Sa <τ ) (Yi (Si ), Yi (t1 ), a, Si , 1)×
i=1
p
Y
fYτ ,Y1 ,Xτ ,T,I(S<τ ) (Yj (τ ), Yj (t1 ), Xj (τ ), τ, 0)
j=1
Z
fY1 ,YS ,XS ,T,I(S≤τ ) (y1 , ys , a, s, 1) = E1 E2 {E3 − C1 E4 } dx1 (5.1)
x1
where
1 1
E1 = fQ1 ,Qs (y1 − cx1 , ys − ca) = −(k+1)/2 −1/2
exp − q c Vf−1 q 0c
(2π) |Vf | 2
2 2 2 2
y1 − cx1 − t1 (νY − cνX ) σY t1 (1 − ρ ) σY t1 (1 − ρ )
qc = Vf =
ys − ca − s(νY − cνX ) σY2 t1 (1 − ρ2 ) σY2 s(1 − ρ2 )
69
(a − x1 − νX (s − t1 )2 )
(a − x1 )
E2 = fS̃a−X (s−t1 ; a−x1 ) = 2
exp − 2
1 (2πσX (s − t1 )2 )1/2 2σX (s − t1 )
(x1 − νX t1 )
E3 = fX1 (x1 ) = 2
(2πσX t1 )−1/2 exp − 2
2σX t1
(2a − x1 − νX t1 )
E4 = fX1 (2a − x1 ) = 2
(2πσX t1 )−1/2 exp − 2
2σX t1
2νX (x1 − a)
C1 = exp 2
σX
Proof. Using arguments from chapter 4, for failed devices we can drop the indicator
variable. The joint density can then be expressed as:
R
= x1
fY1 ,YS ,X1 ,XS ,S (y1 , ys , x1 , a, s) =
R
= x1
fY1 ,YS |X1 ,XS ,S (y1 , ys |x1 , a, s)∗ (5.2)
The first factor in (5.2) is the term E1 , and is equal to fQ1 ,Qs (y1 − cx1 , ys − ca)
by lemma 4.3 and lemma 4.4. The second factor is given by:
The first factor in equation (5.3) represents the term E2 in equation (5.1), and
is given by equation (4.13). The second factor in equation (5.3) is given by equation
(4.21), and represents the term {E3 − C1 E4 } in equation (5.1).
70
5.4.2 Contribution to likelihood from surviving devices
For surviving devices we observe (y1 , yτ , xτ ), where y’s are the observations on
the marker at times t1 and τ respectively, (t1 < τ ), and xτ < a. Time t1 is again
fixed and known. Under the TMDL1 data-structure, the joint density for a survived
device is given by:
fY1 ,Yτ ,Xτ ,I(S>τ ) (y1 , yτ , xτ , 1) = (5.4)
R 2νX (xτ − a)
f (q , q ) fX(τ −t1 ) (xτ − x1 ) − fX(τ −t1 ) (2a − xτ − x1 )exp
x1 Q1 ,Qτ 1 τ 2
×
σX
2νX (x1 − a)
fX1 (x1 ) − fX1 (2a − x1 )exp 2
dx1
σX
(5.5)
where
1 1 −1 0
fQ1 ,Qτ (q1 , qτ ) = exp − q c Vs q c
(2π)−(k+1)/2 |Vs |−1/2 2
q1 = y1 − cx1
qτ = yτ − cxτ
and, k=1, is the number of marker observations, Vs is the variance covariance matrix.
For example, the upper right element of Vs is given by:
2 2 2 2
σY t1 (1 − ρ ) σY t1 (1 − ρ )
Vs =
σY2 t1 (1 − ρ2 ) σY2 τ (1 − ρ2 )
y1 − cx1 − t1 (νY − cνX )
qc =
yτ − cxτ − τ (νY − cνX )
71
(xτ − x1 − νX (τ − t1 )2
2 −1/2
fXτ −t1 (xτ − x1 ) = (2πσX (τ − t1 )) exp − 2
2σX (τ − t1 )
(2a − xτ − x1 − νX (τ − t1 )2
2 −1/2
fXτ −t1 (2a − xτ − x1 ) = (2πσX (τ − t1 )) exp − 2
2σX (τ − t1 )
2
(x 1 − ν t
X 1 )
fX1 (x1 ) = (2πσX 2
t1 )−1/2 exp − 2
2σX t1 )
Proof.
It follows directly from lemma 4.4 that the first factor in equation (5.6) is
equal to fQ1 ,Qτ (q1 , qτ ) ∼ MVN . The second factor is given by:
By theorem 4.2, both factor above have the form of equation (4.21), each given by:
For both failed and survived devices, the joint density between marker obser-
vation time-points should have the same form, as we see in the proof above. The
only contributing factor that should differ (in the likelihood) is the one between
the last marker observation time and the failure-time. This term, as we showed
is a re-started Inverse-Gaussian factor. This repetitive structure will help set the
framework for later deriving the more general longitudinal joint densities. Figure
5.1 illustrates this result.
72
X
x1
t
t1 t1 s
X
x1
xt
t t
t1 t1
Fig. 5.1: Illustration of an observation on the degradation process under a longitudinal
data structure. Rectangular boxes represent density factors given by equation
(4.21), and the oval shape the factor given by equation (4.13)
q
Y
Lθ = fYS ,Y2 ,Y1 ,XS ,T,I(Sa <τ ) (Yi (Si ), Yi (t2 ), Yi (t1 ), a, Si , 1)×
i=1
p
Y
fYτ ,Y2 ,Y1 ,Xτ ,T,I(S<τ ) (Yj (τ ), Yj (t2 ), Yj (t1 ), Xj (τ ), τ, 0)
j=1
73
expressions under the TMD and TMDL1 data-structures.
The derivation of the joint densities in the TMDL1 data-structure call for
integration over latent variable X1 . When more than one intermediate marker is
observed, higher dimensional nested integrals will be needed. High dimensional in-
tegration is not only undesirable from an approximation perspective, it also poses a
serious computational problem. In this section, therefore, we derive a computation-
ally more efficient alternative.
Under the TMDL2 data-structure, the joint density for a failed device is given
by:
fY1 ,Y2 ,YS ,XS ,S (y1 , y2 , ys , a, s) =
R R
= x2 x1 E(s)G1 (x2 )G2 (x1 )H1 (x2 )H2 (x1 , x2 )H3 (x1 )dx1 dx2 = (5.7)
R R
= E(s) x2 G1 (x2 )H1 (x2 ) x1 G2 (x1 )H3 (x1 )H2 (x1 , x2 )dx1 dx2
Definition 5.1. Let the correlation coefficient between Qi and Yj as ρQY , given by:
74
σQs 2 2
E(s) ∼ N µQs + ρQs Y2 (y2 − µY2 ), (1 − ρQs Y2 )σQs
σY2
µQs = E(Ys − cXs ) = s(νY − cνX )
µY2 = νY t2
p
σQs = V ar(Ys − cXs ) =
1/2 1/2
= (V ar(Ys ) + c2 V ar(Xs ) − 2cCov(Ys , Xs )) = (σY2 s(1 − ρ2 ))
√
σ Y2 = σ Y t 2
Cov(Qs , Y2 ) Cov(Ys , Y2 ) − cCov(Xs , Y2 )
ρ Q s Y2 = = =
σQs σY2 σQs σY2
t σ 2 (1 − ρ2 )
= p2 Y √ ⇒ ρQs Y2 = (t2 (1 − ρ2 ))1/2 s−1/2
2
σY s(1 − ρ )σY t2
Then,
where
2
fXt2 −t1 (x2 − x1 ) ∼ N µXt2 −t1 , σX t2 −t1
75
and similarly,
2
H3 (x1 ) = fX1 (x1 ) − exp {2νX (x1 − a)/σX } fX1 (2a − x1 )
2
fX1 (x1 ) ∼ N (νX t1 , σX t1 )
Proof.
where Qs = Ys − cXs = Ys − ca, and c = ρσY /σX . Note that Q(·) is not independent
76
of Y (·). By the same logic, we have:
First factor above is given by H1 (x2 ), the second by H2 (x1 , x2 ) and third by
H3 (x1 ).
Z Z
E(τ ) G1 (x2 )H1 (x2 ) G2 (x1 )H3 (x1 )H2 (x1 , x2 )dx1 dx2 (5.9)
x2 x1
77
The factors in (5.9) are expressed as follows:
σQτ
F (τ ) ∼ N µQτ − ρQτ Y2 (y2 − µY2 ), (1 − ρ2Qτ Y2 )σQ
2
τ
σ Y2
µQτ = E(Yτ − cXτ ) = τ (νY − cνX )
µY2 = νY t2
1/2
σQτ = (σY2 τ (1 − ρ2 ))
√
σY2 = σY t2
Cov(Qτ , Y2 ) cov(Yτ , Y2 ) − cCov(Xτ , Y2 )
ρQτ Y2 = = =
σQτ σY2 σQτ σY2
t σ 2 (1 − ρ2 )
= p2 Y √
σY τ (1 − ρ2 )σY t2
⇒ ρQτ Y2 = (t2 (1 − ρ2 ))1/2 τ −1/2
and,
Proof.
78
Due to theorem 4.2 and lemma 4.3 we get:
79
In a general case, under the GENL data-structure, for a failed device, the joint
density is given by:
fY1 ,...,YK ,YS ,XS ,T,I(S<τ ),K (y1 , . . . , yk , ys , a, s, 1, k) =
R R
xk
... fQ1 ,...,Qk ,Qs (q1 , . . . , qk , qs )fS (s − tk ; a − xk )∗
x1
k
Y 2νX (xj − a)
∗ fXj −Xj−1 (xj − xj−1 ) − exp 2
fXj −Xj−1 (2a − xj − xj−1 )
j=1
σX
(5.10)
Proof.
(a) = fQ1 ,...,QK |X1 ,...,XK ,XS ,S,K (q1 , . . . , qk |x1 , . . . , xk , a, s, k) (5.12)
Let X k−1 = {X1 , . . . , Xk−1 }, then factor (b) in (5.11) is written as:
80
The randomness of K = k is entirely captured by S = s. Therefore
(b1) = fS̃a−x (s − tk ; a − xk )
k
Take one term in factor (b2). Due to theorem 4.2 we only condition on the last
observation on X at time tj−1 . Then,
Therefore, putting things together, and plugging in factor (b1) and (b2) into
(b) we get:
(b) = fS̃a−x (s − tk ; a − xk )∗
k
Qk 2νX (xj − a)
∗ j=1 {fXj −Xj−1 (xj − xj−1 ) − exp 2
fXj −Xj−1 (2a − xj − xj−1 )}
σX
where x0 = 0, t0 = 0.
81
fixed and non-random. In the general case, under the GENL data-structure, the
joint density for surviving device is given by:
Proof.
fY1 ,...,Yn ,Xn ,I(S≥τ ) (y1 , . . . , yn , xn , 1) =
R R
= xn−1 . . . x1 fY n |X n ,I(S≥τ ) (y|x, 1) ∗ (a)
∗ fX n ,I(S≥τ ) (x, 1)dx1 . . . dxn−1 (b)
Factor (b) is expanded like before, using theorem 4.2 we can proceed as follows:
n
Y 2νX (xj − a)
(b) = fXj −Xj−1 (xj − xj−1 ) − fXj −Xj−1 (2a − xj − xj−1 )exp 2
j=1
σX
82
5.7 Summary and Conclusions
We were able to identify the two density factors needed to analytically express
the joint densities for both failed and surviving devices in the likelihood function.
These are namely, (1) the re-started inverse-Gaussian fS̃a−X(t ) (s − tj ; a − xj ), given
j
by equation (4.13) and (2) fXj ,I(S≤tj ) (xj , 0), given by equation (4.21) .
Although the analytical structure is simple, larger multidimensional marker
observations require an equal number of integrations, a fact that will burden com-
putations. Some simplification can be achieved, as we show, by factoring out terms
that are not functions of the space over which we integrate. However, high dimen-
sional nested integrals still remains a computational limitation.
83
6. CASE STUDIES - TERMINAL AND LONGITUDINAL
Table 6.1 presents the MLE notations under the four data-structures we ex-
amine in this chapter.
• At what sample sizes are estimators under the two-data-structures finite sam-
ple efficient? In other words, at what sample size does the empirical variance
approach the Cramer Rao lower bound? Is this sample size different for the
two data-structures?
• What is the Asymptotic Relative Efficiency (ARE) between the two data-
structures? In other words, for a large enough sample size, what is the ratio
of the variances for the corresponding data-structures.
85
• Even though we do not expect MLEs to be efficient under small samples, the
relative performance of the two estimators for small sample sizes is of interest.
We calculate therefore, and analyze the ARE for small samples
definition 3.9 is defined as the ratio of estimated asymptotic variances: V̂ (µ̂2 )/V̂ (µ̂1 ).
The asymptotic variance of µ̂, V (µ̂), is derived through the Delta method, which
√ D
→ N (0, σ 2 ), where νX and σ 2 are finite valued constants
says that if n(νX − ν̂X ) −
D
and −
→ denotes convergence in distribution, then
√ D
→ N (0, g 0 (νX )2 σ 2 )
n(g(νX ) − g(ν̂X )) −
V (g(ν̂X )) = V (µ̂) = g 0 (νX )2 σ 2
In the ratio defined by the ARE, the term g 0 (νX )2 cancels out, and we are left with
86
the ratio of asymptotic variances for νX of TMD vs TM.
For the MLE’s obtained under each data-structure, the estimated asymptotic
−1
variance nV̂ (νˆX ) is the upper-left element (Iobs )1,1 of the inverse of the observed
information matrix [Remarks 3.2 & 3.3], averaged over R simulations to obtain the
r
estimated asymptotic relative average efficiency ARE
\R , with ν̂X the MLE for νX ,
obtained under the rth simulated data-set.
1 PR TM
r=1 V̂ (ν̂X )r
[ X , ν̂X ) = R
ARE(ν̂ T M D T M
(6.1)
1 PR T MD
V̂ (ν̂X )r
R r=1
87
ing phenomena: degradation of the device causes changes in the marker variable
distribution. Our inference model, however, works in the other direction: it takes
marker observations and infers the distribution of the degradation variable at each
time point. Next we describe the simulation model.
Physically the rationale behind the selected parameter values can be motivated
from the degradation of solder joints. Solder joints hold electronic components (like
a capacitor, or a resistor) to a printed circuit board, and can, with use (abuse),
develop micro-cracks, which in time can grow and cause reliability problems. These
cracks can grow or shrink depending on usage or environmental conditions. Then,
X(t) in the model can be thought of as the length or size of the crack. It is not
unreasonable to have high variance in this process, especially since the crack can
entirely close ”heal” (given the right conditions), and then snap back to a fully
”opened” state. So the signal to noise ratio (νX /σX ) is reasonably less than 1.
The drift and variability of the marker process are not as important (at least in its
physical interpretation), but the correlation coefficient ρ of course is. We see the
influence of varying ρ in studying the ARE of TMDL1 vs TMD. We also consider
other parameter combinations to study more general patterns in estimation.
Before getting to the simulation results, we first define the simulation design.
We discuss the method of generating the degradation, marker and lifetime data
for the TM data-structure. We start by partitioning the total time-on-test [0, τ ]
into N = 200 equally sized time intervals ∆t, and construct a time vector ts =
(ts1 , . . . , tsN ), such that s s
P
∀i ∆t = tN = τ , i = 1, . . . , N . We then generate an
88
with common parameters, as:
∆X(ti ) ∼ N νX ∆ts , σX
2
∆ts
(6.2)
89
from Iobs compared to the empirical sample variance, under both the TM and
TMD data-structures. The close agreement between the Iobs -derived and empiri-
cal columns, especially for large samples, is predicted by asymptotic MLE theory.
As expected the uncertainty in the drift parameter estimates are consistently smaller
under the TMD data-structure. In particular the large-sample Cramer-Rao lower
bound is smaller for the TMD than it is for TM. Both TM and TMD seem to attain
efficiency, up to the accuracy of the simulation, at around a sample size of 100 (This
can be seen by comparing Cramer Rao lower bound to the empirical variances).
Tab. 6.2: Asymptotic Standard Errors for νX under TM and TMD, from Observed Infor-
mation vs. Empirical
Obtained from Iobs Empirical
Sample Size (n) TM TMD TM TMD
20 0.01495 0.01207 0.01604 0.0127
40 0.01052 0.00866 0.01053 0.00889
80 0.00750 0.00616 0.0077 0.00616
160 0.00530 0.00437 0.00554 0.00446
320 0.00377 0.00309 0.00383 0.00317
640 0.00266 0.00219 0.00262 0.00219
1280 0.00188 0.00155 0.0019 0.00157
2560 0.00133 0.00110 0.00138 0.00108
Table 6.3 tabulates the central limit theorem based confidence intervals (CI)
for the empirical asymptotic variances in table 6.2, and provides a measure of how
precise the agreement is between Iobs -derived and empirically-derived columns in
r
table 6.2. MLE’s ν̂X , r = 1, . . . , R, are independent and asymptotically normal, or
we can say approximately normal for large enough n. For large enough n, therefore,
r
we have that ν̂X ∼ N (νX , V (νX )), where V (νX ) is the asymptotic variance of νX .
− ν̂¯X )2 /(n − 1),
P r
The estimated asymptotic variance is given by: V̂ (νX ) = r (ν̂X
where ν̂¯X is the empirical average of ν̂X over R samples, and (n − 1)V̂ (νX )/V (νX ) ∼
χ2n−1 . Therefore, a 100(1 − α)% CI for V (νX ) is given by:
n−1 n−1
2
V̂ (νX ) < V (νX ) < 2 V̂ (νX )
χα/2,n−1 χ1−α/2,n−1
90
We are interested in the CI for the standard deviations (SD), so we take the
square root of the upper and lower limit of the CIs for the variance. Therefore we
have: s s
n−1 d n−1
SD(νX ) < SD(νX ) < SD(ν
d X) (6.3)
χ2α/2,n−1 χ21−α/2,n−1
q
where SD(ν
d X) = V̂ (νX )
Tab. 6.3: Central Limit Theorem-based confidence intervals for asymptotic empirical Stan-
dard Deviations under TM and TMD
TM TMD
Sample
Size (n) Upper Lower Upper Lower
20 0.01219 0.02342 0.00965 0.0185
40 0.00862 0.01352 0.00728 0.01141
80 0.00666 0.00912 0.00533 0.00729
160 0.00499 0.00622 0.00401 0.00501
320 0.00355 0.00415 0.00294 0.00343
640 0.00248 0.00277 0.00207 0.00231
1280 0.00182 0.00197 0.00151 0.00163
2560 0.00134 0.00141 0.00105 0.00111
In failure tests, engineers are often interested in reducing test durations and
accelerating factors/conditions. Longer test durations cost more money and re-
sources, and overly accelerating conditions can change targeted failure mechanisms.
Failure test-designs that use low or no accelerating conditions are preferred. In
addition, correlation between degradation X and the marker Y is often weaker
than expected. There is interest, therefore, to investigate the efficacy of the TMD
data-structure when τ is small and correlation ρ is weak. A simulation experi-
T MD TM
ment computes ARE(ν̂
[ X , ν̂X ) for pairs of (ρ, τ ) from ρ = (0, 0.3, 0.6, 0.9) and
τ = (5, 10, 15, 20). The correlation ρ is varied from weak to strong and τ from short
to long.
For each pair (ρ, τ ), and fixed n = 500, we simulated R = 100 data-sets and
∗
computed MLEs θ̂ r = (νX , νY , σX , σY )r , and their asymptotic variance estimates
91
Tab. 6.4: Asymptotic relative efficiency of µ̂T M D versus µ̂T M for different (ρ, τ ) combina-
tions, with θ = (0.1, 1.0, 0.2, 0.1)
ρ/τ 4 7 10 13
0 5.9614 2.1593 1.4759 1.2488
0.3 5.6893 2.1465 1.4505 1.2336
0.6 5.4270 2.0464 1.4106 1.2088
0.9 4.0123 1.6320 1.2355 1.1262
∗
V̂ (θ̂ r ), r = 1, . . . , R, under both TM and TMD data separately. We calculate
ARE(·),
[ according to equation (6.1) by averaging V̂ (·)r . Table 6.4 shows the ARE
results from an experiment with θ ∗ = (0.1, 1.0, 0.2, 0.1). In table 6.4 we point out
the improvement in inference for low (ρ, τ ) combinations. This result indicates
that under the bivariate Wiener model the TMD data-structure improves inference
under smaller sample sizes and weaker correlation coefficients. This result is very
promising because it gives preliminary justification for a TMD data-structure.
Tab. 6.5: Asymptotic relative efficiency of µ̂T M D versus µ̂T M for different (ρ, τ ) combina-
tions, with θ = (0.1, 1.0, 0.4, 0.1)
ρ/τ 4 7 10 13
0 2.44236 1.68550 1.43906 1.31727
0.3 2.40830 1.67795 1.41769 1.31341
0.6 2.29307 1.59947 1.394463 1.29086
0.9 1.80918 1.40500 1.263771 1.19975
and with a different θ ∗ = (0.1, 1.0, 0.4, 0.1). We can see that the pattern of AREs
across the (ρ, τ ) grid is qualitatively similar. However, because the variance of the
92
simulated Wiener process is larger, now σX = 0.4, it causes a larger proportion of
devices to fail, which increases the ARE for each (ρ, τ ) combination.
The above results in tables 6.4 and 6.5 can be explained as follows. With
small τ ’s there are more censored lifetime observations in each simulated dataset,
which means there are more observations on the terminal degradation in the TMD
data-structure. With more observations on degradation, in turn, we expect im-
proved inference under the TMD data-structure, and therefore higher AREs. AREs
greater than 1 indicate the efficacy of the TMD data-structure, and AREs close to
1 indicate that the two data-structures provide the same inference power. When
θ ∗ = (0.1, 1, 0.2, 0.1), the percentage of failed devices is: 11%, 34%, 74%, and 97%,
respectively for the four τ levels. When θ ∗ = (0.1, 1.0, 0.4, 0.1) the percentage of
failed devices is: 33%, 68%, 89%, and 99%, respectively for the four levels of τ . We
see that with higher variance in the process, we have more failures under each τ , and
therefore, access to less degradation information on surviving device, and therefore,
as expected, lower AREs.
No. νX νY σX σY No. νX νY σX σY
1 0.1 0.9 0.2 0.1 12 0.1 1.0 0.2 0.4
2 0.1 0.5 0.2 0.1 13 0.1 1.0 0.2 0.8
3 0.1 0.1 0.1 0.1 14 0.1 1.0 0.2 1.0
4 0.1 2.0 0.2 0.1 15 0.1 1.0 0.2 1.5
5 0.1 4.0 0.2 0.1 16 0.1 1.0 0.2 2.0
8 0.1 1.0 0.2 0.1 17 0.4 1.0 0.2 0.1
6 0.1 1.0 0.4 0.1 18 0.8 1.0 0.2 0.1
7 0.1 1.0 0.8 0.1 19 1.0 1.0 0.2 0.1
9 0.1 1.0 1.0 0.1 20 1.5 1.0 0.2 0.1
10 0.1 1.0 1.5 0.1 21 2.0 1.0 0.2 0.1
11 0.1 1.0 2.0 0.1
More generally, we compute the AREs for the same set of (ρ, τ ) pairs, failure
threshold level a = 1, and time increment ∆t = 0.01, using various combinations of
93
2.8 2.8
2.6 2.6
2.4 2.4
2.2 2.2
2 sX 2 sX
ARE
ARE
1.8 1.8
1.6 1.6
1.4 1.4
1.2 1.2
1 1
4 7 10 13 4 7 10 13
Fig. 6.1: General ARE patterns of µ̂T M D versus µ̂T M , under increasing σX
94
7 50
r=0 45 r = 0.6
6
40
35
5
30
sY
ARE
ARE
4 25
sY
20
3
15
2 10
1 0
4 7 10 13 4 7 10 13
t t
Fig. 6.2: General ARE patterns of µ̂T M D versus µ̂T M , sunder increasing σY
for ρ = 0 (left) and ρ = 0.6 (right). Here we see that under fixed τ , νY does not
influence estimation. The results from figures 6.2 and 6.3, indicate that inference
is not affected by the marker magnitude at termination, but rather its variance.
Although both data-structures include terminal marker observations, inference suf-
fers more under TM than TMD because, under the TM data-structure, inference
is only based on the marker observations, whereas under the TMD data-structure,
inference is based on both marker and degradation data.
Increasing the drift of the degradation process on the other hand has more
prominent consequences on AREs. Figure 6.4 plots the general ARE patterns for
ρ = 0 (left) and ρ = 0.6 (right), under increasing νX , as a function of τ . As expected,
under each fixed τ , larger νX , increases the ARE. This is highlighted especially for
small values of τ . These results make sense, because under high drifts, we expect
a larger proportion of failed devices, and therefore as per our hypothesis, improved
inference under the TMD as opposed to the TM data-structure.
Generally, these results validate our hypothesis that inference improves with
enhanced data, and that the TMD data-structure is more efficient at ML-estimation
under small failure-time samples. These results can be used to justify reducing
95
7 6
4.5
5
4
ARE
ARE
3.5
nY nY
3
3
2.5
2
2
1.5
1 1
4 7 10 13 4 7 10 13
Fig. 6.3: General ARE patterns of µ̂T M D versus µ̂T M , under increasing νY
6 5.5
r=0 5 r = 0.6
5
4.5
4
4
3.5
nX nX
ARE
ARE
3 3
2.5
2
2
1.5
1
1
0 0.5
4 7 10 13 4 7 10 13
Fig. 6.4: General ARE patterns of µ̂T M D versus µ̂T M , under increasing νY
96
accelerating conditions in accelerated failure tests. By reducing the accelerating
conditions, and keeping the test-duration τ fixed, we are likely to observe fewer
failed devices by time τ . However, based on the above results, we need a smaller
failure-time sample in the TMD data-structure to get the same predictive power
as the TM data-structure, and therefore we can ”afford” to reduce the designed
accelerating conditions in the planned failure-test. From an engineering perspective,
reducing accelerating conditions in failure tests is a welcomed option, because, we
we discuss in the introduction, highly accelerated test conditions may alter targeted
failure-mechanisms.
These results can also be used to justify shorter test-durations for planned
failure-tests. Typically, failure-tests are conducted over a time-period long enough
to see enough devices fail. Shorter test-durations are not only less costly, they are
sometimes required due to short product life-cycles.
High values of ρ make the marker data more closely associated to the failure
mechanism, and therefore we expect lower AREs when ρ approaches 1. When
the marker is perfectly correlated to the degradation variable, any observations on
the degradation should not enhance inference. This argument is validated in the
simulation [Link] general we notice that for stronger correlations we get smaller
ARE. For long test-times, most devices fail, and the TMD data-structure, therefore,
looses valuable degradation information from surviving devices.
The improvement in estimation under the TMD data-structure is due to obser-
vations on terminal degradation. In tests where degradation is latent, as motivated
in this thesis, we argued and proved with the results above that access to termi-
nal degradation data reduces estimation variance in a bivariate Wiener model. In
failure-test laboratories, terminal degradation is often measured using expensive
equipment or proprietary techniques, and is therefore valuable information. Un-
der latent degradation conditions, the above results show the efficacy of terminal
97
degradation data in reducing estimation variance. These results therefore, suggest
investing in failure-analysis equipment capable of measuring terminal degradation.
Predictions are computed using the predictive inference equations (4.35) and
(4.36) and the MLEs as plug-in estimators. Given survival and a marker obser-
vation at a sequence of time points ti , Figure 6.5 plots the predicted conditional
degradation (left) and future failure-time densities (right). The probability density
function of the degradation variable at time ti is plotted by connecting probabilities
computed for a range of degradation values. Similarly the probability density func-
tion of the failure-time variable at time ti is plotted by connecting the probabilities
computed for a range of future failure-times. Predictions qualitatively capture the
behavior/trend of the latent degradation process, and the decreasing uncertainty in
failure-time predictions.
The shape of the density is plotted by connecting point probabilities evaluated
for a vector-valued dependent variable. For example, we vary the level of degrada-
tion from -1 to 1 in small increments and evaluate the probability specified by the
predictive inference equations for each xi
98
Fig. 6.5: Predicted degradation and survival density as a function of time
For each sample size n, we simulate sample observations with the parameter
set θ = (νX , νY , σX , σY , ρ) = (0.1, 1.0, 0.4, 0.1, 0.75). Lifetimes were generated in
the same way as before (see section 6.1), with ∆t=0.01, a=1, and τ =10. In each
simulation run j, the intermediate marker observation is made at time-point t1 =
τ /2. This way the marker observation time is independent of the failure-time and
always less than the end-of-test time. When t1 > s the contribution of that failed
device to the likelihood is made through the TMD data-structure.
Table 6.7 reports the asymptotic standard errors for νX computed from Iobs
compared to the empirical sample variance-covariance matrix, under both the TMD
and TMDL1 data-structures. We observe again a close agreement between the
Iobs -derived and empirical columns. We observe, contrary to our expectation no im-
provement in estimation results under the TMDL1 data-structure in comparison to
the TMD data-structure. This result is likely a result of simulation errors. However,
under TMDL2 we actually see improvement in estimation as we discuss in the next
section.
Asymptotic relative efficiency is used to test the efficacy of TMD vs TMDL1.
99
Tab. 6.7: Asymptotic Standard Errors for νX under TMD and TMDL1, from Observed
Information vs. Empirical
Obtained from Iobs Empirical
Sample Size (n) TMD TMDL1 TMD TMDL1
20 0.00711 0.00621 0.00729 0.00828
40 0.00511 0.00533 0.00536 0.0058
80 0.00366 0.00369 0.00354 0.00354
160 0.00261 0.00263 0.00272 0.00272
320 0.00185 0.00187 0.00201 0.00201
640 0.00131 0.00131 0.00127 0.00127
Tab. 6.8: Central Limit Theorem-based confidence intervals for asymptotic empirical vari-
ances under TMD and TMDL1
TMD TMDL1
Sample size (n) Lower Upper Lower Upper
20 0.00554 0.01064 0.00629 0.01209
40 0.00439 0.00688 0.00475 0.00744
80 0.00306 0.00419 0.00306 0.00419
160 0.00245 0.00305 0.00245 0.00305
320 0.00186 0.00217 0.00186 0.00217
640 0.00120 0.00134 0.00120 0.00134
100
dation. In this case it matters how strongly correlated the marker variable is to the
degradation, and table 6.9 verifies this hypothesis.
Once again, we have grounds to advertise the TMDL data-structures in ap-
plication areas with few failure observations. Although direct information on the
degradation variable is shown here to ”outweigh” information on the marker vari-
able (under certain conditions), it remains to be seen if there comes a point when
this is no longer the case. In other words, for estimation purposes, does access to
terminal degradation information become irrelevant when we have a lot of longitu-
dinal marker information? One must also consider computational complexity and
computational efficiency to make the comparison fair, after all, multivariate longitu-
dinal marker observations require more complex multivariate (nested) integration.
Tab. 6.9: Relative efficiency of µ̂T M DL1 versus µ̂T M D , for different (ρ, τ ) combinations
ρ/τ 4 7 10 13
0 1.00053 0.99984 0.99974 0.99983
0.3 1.01364 1.00681 1.00269 0.99854
0.6 1.07101 1.04112 1.01713 0.99849
0.9 1.24527 1.16814 1.06883 1.01328
In effort to gain further insight into the contribution that additional marker
observations have on estimation, we investigate the TMDL2 data-structure. We
again simulate R=250 independent data-sets D, for n = (20, 40, 80, 160), and
computed the MLEs θ̂ for each. MLEs from the TMD data-structure are repre-
sented by vector θ̂ T M D , from the TMDL1 data-structure by θ̂ T M DL1 , and from the
TMDL2 data-structure by θ̂ T M DL2 . Under TMDL2, for failed devices, we observe
Y = (Y (t1 ), Y (t2 ), Y (s)), and for surviving devices Y = (Y (t1 ), Y (t2 ), Y (τ )), and
as always 0 ≤ s ≤ τ , and t1 < t2 < τ . Devices that fail before t2 contribute to the
101
likelihood function using the TMDL1 data-structure, and devices that fail before t1
contribute through the TMD data-structure.
The same parameter set θ = (0.1, 1.0, 0.4, 0.1, 0.75) is used, and lifetimes are
generated in the same way as before, with ∆t=0.01 and a=1. In each simulation
run j, the intermediate marker observations are made at t1 = τ /3 and t2 = 2τ /3.
Table 6.10 reports the asymptotic standard errors for νX computed from Iobs . The
TMDL2 data-structure, on average, results in smaller asymptotic standard errors,
than the TMDL and TMD structures, indicating therefore, without the aid of AREs,
improved estimation. Further higher dimensional longitudinal data-structures are
not investigated in this work, due to lack of computing power for evaluating higher
dimensional integrals. Using the results in table 6.10 however, we can assume that
inference will improve with more marker observations. It is not clear from these
results if the influence of additional marker observations will attenuate.
Tab. 6.10: Asymptotic Standard Errors for νX under TMD, TMDL1 and TMDL2 from
Observed Information
Obtained from Iobs
Sample Size TMD TMDL1 TMDL2
20 0.00612 0.00612 0.00604
40 0.00483 0.00483 0.00466
80 0.00365 0.00365 0.00360
160 0.00270 0.00270 0.00193
102
We chose ST 3 − 0 arbitrarily amongst the 6 settings.
One of the key features in the data is a latent degradation variable, which
reflects the fact that the degradation variable is not observed. Because each device is
tested to failure, however, we can assume that for each failed device, the degradation
variable reached a fixed known threshold at the failure time. This assumption allows
us to apply the FHT models presented in this thesis. In the simplest case, we
simulate degradation to be linear in cycles, with heterogeneous drift and variance
across devices. In the analysis we make the following assumptions:
• Degradation paths are assumed linear and proportional to the device’s lifetime
• If the device is not observed on the 200th cycle, then collect data on the cycle
closest to cycle 200
103
Tab. 6.11: Illustration of ST3-0 data-structure with augmented degradation data
Usage Setting ST3-0
Device Cycles z1 z2 . . . z21 Deg
1 4 491 601 . . . 0.22 0.002
1 6 491 607 . . . 0.19 0.03
1 23 490 607 . . . 0.24 0.009
.. .. .. .. .. .. ..
. . . . . . .
1 31 491 607 . . . 0.24 0.012
1 45 490 608 . . . 0.24 0.018
1 75 496 607 . . . 0.29 0.023
To fit this data into the TM and TMD data-structures, we are interested in
suitable marker variables that can be used. As a first data-cleaning step, variables z1 ,
z5 , z18 and z19 with zero standard deviation are removed from the data-set, leaving
17 covariates. Figure 6.6 plots a colormap surface of the empirical correlation matrix
between each of the 17 covariates and the degradation variable. From examining the
correlation matrix we see that variables 9,13, and 14 are strongly correlated with
each other and therefore redundant. From the remaining variables, only four show
considerable correlation to the degradation variable. We continue our analysis using
covariates z4 , z8 , z11 and z15 .
Using principal component analysis (PCA), we reduce the dimensionality fur-
ther, and ultimately derive one variable that can be used as the marker variable.
PCA forms a new set of uncorrelated (not necessarily independent) variables that we
denote by z 0 . The PCA scores on the most dominant eigenvector show the strongest
correlation to the degradation variable (∼ 0.42), and therefore this variable is used
as the marker variable.
Table 6.12 compares the asymptotic standard errors for TM vs TMD under
increasing sample sizes n = (20, 40, . . . , 218). In these results we observe and val-
idate that under the TMD data-structure, parameters are consistently estimated
with lower uncertainty. Also its interesting to observe that under both models the
estimates for ρ are in close agreement to the empirical correlation between marker
104
s2 s3 s4 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s20 s21 deg
s2 1 0.547 0.643 0.45 0.554 0.515 0.273 0.335 0.674 0.578 0.517 0.18 0.642 0.547 0.557 0.484 0.498 0.612
s3 0.547 1 0.619 0.433 0.537 0.5 0.264 0.323 0.649 0.558 0.495 0.176 0.62 0.509 0.533 0.483 0.481 0.583
s4 0.643 0.619 1 0.496 0.647 0.528 0.235 0.374 0.775 0.676 0.528 0.127 0.729 0.617 0.636 0.565 0.569 0.68
s6 0.45 0.433 0.496 1 0.43 0.373 0.18 0.231 0.525 0.456 0.374 0.109 0.492 0.344 0.428 0.38 0.381 0.502
s7 0.554 0.537 0.647 0.43 1 0.434 0.172 0.322 0.673 0.575 0.435 0.076 0.633 0.545 0.545 0.497 0.495 0.586
s8 0.515 0.5 0.528 0.373 0.434 1 0.87 0.369 0.534 0.436 0.945 0.818 0.569 0.49 0.515 0.411 0.412 0.653
s9 0.273 0.264 0.235 0.18 0.172 0.87 1 0.244 0.217 0.161 0.87 0.948 0.29 0.261 0.275 0.189 0.188 0.44
s10 0.335 0.323 0.374 0.231 0.322 0.369 0.244 1 0.389 0.335 0.371 0.191 0.387 0.338 0.346 0.298 0.284 0.381
s11 0.674 0.649 0.775 0.525 0.673 0.534 0.217 0.389 1 0.716 0.533 0.101 0.764 0.654 0.671 0.593 0.6 0.706
s12 0.578 0.558 0.676 0.456 0.575 0.436 0.161 0.335 0.716 1 0.438 0.057 0.666 0.562 0.576 0.511 0.519 0.613
s13 0.517 0.495 0.528 0.374 0.435 0.945 0.87 0.371 0.533 0.438 1 0.819 0.571 0.494 0.514 0.414 0.41 0.656
s14 0.18 0.176 0.127 0.109 0.076 0.818 0.948 0.191 0.101 0.057 0.819 1 0.183 0.17 0.185 0.102 0.102 0.348
s15 0.642 0.62 0.729 0.492 0.633 0.569 0.29 0.387 0.764 0.666 0.571 0.183 1 0.622 0.635 0.554 0.563 0.688
s16 0.547 0.509 0.617 0.344 0.545 0.49 0.261 0.338 0.654 0.562 0.494 0.17 0.622 1 0.544 0.479 0.481 0.604
s17 0.557 0.533 0.636 0.428 0.545 0.515 0.275 0.346 0.671 0.576 0.514 0.185 0.635 0.544 1 0.497 0.498 0.608
s20 0.484 0.483 0.565 0.38 0.497 0.411 0.189 0.298 0.593 0.511 0.414 0.102 0.554 0.479 0.497 1 0.433 0.533
s21 0.498 0.481 0.569 0.381 0.495 0.412 0.188 0.284 0.6 0.519 0.41 0.102 0.563 0.481 0.498 0.433 1 0.53
deg 0.612 0.583 0.68 0.502 0.586 0.653 0.44 0.381 0.706 0.613 0.656 0.348 0.688 0.604 0.608 0.533 0.53 1
Fig. 6.6: Correlation matrix between covariates (1 through 17) and degradation (18)
and degradation (∼ 0.42 ). Similarly we observe that the MLE for νX is in close
agreement with the empirical rate, which is calculated by 1/(average number of cy-
cles to failure) ∼ 200 = 0.005. Table 6.13 compares the asymptotic standard errors
and MLEs for ρ across the same range of sample sizes. Here again we observe an
improved inference under the TMD data-structure. For larger n, although the vari-
ance is reduced (under both data-structures) the MLE for ρ seems to be relatively
biased to the empirical correlation.
The proportion of failed to survived devices in each sample size is random and
approximately equal to 50% as noted in Table 6.11. This is a result of constructing
a sample by adding records from subsequently observed censored and failed devices,
starting from the 1st device. For example, for a sample of size n = 20, we construct
the data set by adding an additional row to the data-set for each device number,
starting from device 1. In this way, the proportion of devices out of n = 20 that fail
105
Tab. 6.12: Asymptotic Standard Error for νX under TM and TMD from Observed Infor-
mation. Empirical drift ∼ 0.005
MLE SE
Sample Size TM TMD TM TMD Prop. Failed
20 0.00511 0.00526 0.0003445 0.000244 0.55
40 0.005 0.00505 0.0002116 0.0001598 0.525
60 0.00483 0.00497 0.0002000 0.0001305 0.466
79 0.00486 0.00501 0.0001757 0.0001150 0.482
99 0.00483 0.00499 0.000153 9.91E-05 0.484
119 0.00481 0.00498 0.0001385 8.79E-05 0.478
139 0.00485 0.00502 0.0001279 8.18E-05 0.489
159 0.00483 0.005 0.0001228 7.81E-05 0.484
179 0.00482 0.00498 0.0001128 7.22E-05 0.48
199 0.00481 0.00499 0.0001099 6.94E-05 0.482
218 0.0048 0.00499 0.000108 6.75E-05 0.481
or survive is random.
Tab. 6.13: Asymptotic Standard Error for ρ under TM and TMD from Observed Infor-
mation. Empirical correlation = 0.42
MLE SE
Sample Size TM TMD TM TMD
20 0.42394 0.44558 0.18111 0.16280
40 0.39581 0.44455 0.13869 0.11626
60 0.40546 0.41031 0.11424 0.09983
79 0.43247 0.42076 0.09574 0.08643
99 0.43147 0.41128 0.08417 0.07767
119 0.43723 0.40441 0.07613 0.07194
139 0.41361 0.38366 0.07177 0.06783
159 0.40293 0.38431 0.06872 0.06307
179 0.40918 0.39213 0.06449 0.05918
199 0.36231 0.35231 0.06408 0.058395
218 0.35679 0.33982 0.06128 0.056277
106
Fig. 6.7: Relative efficiency of µ̂ from G2 vs. G1 for different (ρ, n) combinations evaluated
with gas-turbine engine degradation data
Figure 6.4 plots the results of the ARE experiment, and tells a different
story from the ARE results in the simulations. In this case, the proportion of
failed to surviving devices remains the same for all sample sizes. We observe that
the efficacy of the TMD data-structure decreases with increasing ρ and decreasing
sample size n. The rationale behind these results can be explained in terms of
information content. When ρ is high it favors the TM data-structure because it can
draw more information on surviving devices from the marker observations. From
the TMD perspective, when ρ is high, therefore the marker information becomes
more valuable, it takes away the advantage of observing terminal degradation. On
the other hand, with large sample sizes the TMD data-structure slowly regains its
advantage.
107
7. COVARIATES AND REGRESSION STRUCTURES
Emphasis in this thesis has been given to the development of parametric infer-
ence procedures for the bivariate Wiener model under various extensions to the TM
data-structure. The effect of covariates still remains to be incorporated into these
models. Covariates are measurable variables that contain information about device
specific performance and environment. Covariate data in an engineering setting
are usually collected from sensors embedded onto or near the device, and designed
to measure targeted information related to the health/degradation of the device.
Reliability models that incorporate covariate data together with lifetime data are
anticipated to improve reliability predictions.
Much literature has been devoted to reliability models for lifetime data with
covariates. This is especially true in what are known as accelerated failure time
models, with early work by Epstein and Sobel 1963, Singpurwalla 1970, and many
others, some of which we discuss later in this chapter. The main purpose of this
chapter is to discuss the relevance of using covariates in degradation models, and
specifically in FHT models.
Lee and Whitmore 2006 introduce Threshold Regression models, which are
FHT models that include regression structures. Regression structures allow effects
of covariates to explain some or all of the dispersion of the data, thereby taking
account of variability ad sharpening inferences [Whit, Lee 2006]. The idea of using
covariates to aid estimation and inference can be found in a broad range of literature
in survival and reliability analysis. The main idea is that unknown distribution or
process parameters are reparameterized as functions of covariates.
In FHT models, one can think about reparameterizing the drift parameter νX
as νX = α − βZ, where Z is a vector valued covariate, and α, β are new unknown
parameters that need to be estimated. In electronics reliability experiments, ac-
celerated lifetimes are most common, induced by higher than normal experimental
stress conditions, such as high temperature, humidity and pressure. It becomes
important to model the influence of covariates that measure the stresses and other
environmental conditions, on the rate of degradation. In other words, they include
the effect that covariate information has on the drift of the degradation process.
This relationship is especially important during predictive inference calculations on
test devices. In the example used above, instead of using the MLE νˆX as a plug-in
estimator in predictive inference equations, we can use α̂ − β̂Zi for device i, there-
fore explaining the drift accounting for heterogeneous covariate observations across
the training sample.
One question that arises is what regression structure to use that best captures
the functional relationship between covariates and the dependent variable of interest,
like, drift or more generally, degradation. In cases where failure-mechanisms are well
understood, then PoF models can motivate the appropriate regression structure.
PoF models typically relate covariates at a point in time to either the degradation
level or to the lifetime scale parameter. In case where lifetimes are attained under
what are called accelerated test conditions, PoF models also try to account for
the acceleration. When failure-mechanisms are not well understood, then data-
driven approaches can help identify the statistical relationship. In the remainder of
the chapter we discuss various regression structures that can be considered in the
context of FHT models.
109
7.1 Multiplicative Hazards Regression Models
In survival analysis the hazard rate λ(t) captures the instantaneous probability
of failure at time t, and forms the basis for multiplicative hazard regression models:
where λ0 (t) is a baseline hazard function and r(z i (t)) is a positive valued device-
specific relative-risk multiplier. Typically, the relative-risk is specified by r(z i (t)) =
exp(θ T z i (t)) with parameter vector θ. The hazard ratio between any two devices,
λi (t)/λj (t) = exp(θ T (z i (t) − z j (t))), is constant if the difference in covariates is
constant in time, specifying therefore a proportional hazard model.
In Cox’s proportional hazard regression model λi (t) = λ0 (t)exp(θ T z i (t)) ,
the baseline hazard is left unspecified. Its semi-parametric form can help reduce
estimation bias of covariate effects. The parameter vector θ can be estimated using
the partial likelihood [74].
Accelerated tests are typically used to collect information on the life distribu-
tion or performance over time of products. Meeker and Escobar 1993 provide an
excellent review of research and issues in accelerated testing. In their 2002 book,
Bagdonavicius and Nikulin present a comprehensive review of univariate accelerated
life models, Nelson 1990 describes accelerated life models and life-stress relationships
such as the Arrhenius, the inverse-power, and the fatigue relationships.
Accelerated failure time regression models, provide an alternative to the com-
monly used proportional hazards models. AFT regression models assume that the
effect of a covariate is on the failure-time itself. AFT regression models are typi-
cally based on log transformations of the failure-time. For example, in a log-linear
110
failure-time model log(si ) = µi + β T z i + , with −∞ < µi < ∞, β > 0 and is a
normally i.i.d. error variable. The survival function and hazard rate are then given
by:
S(t|z i ) = S0 (texp(−β T z i ))
In FHT models, marker variables, as discussed earlier, form the basis for in-
ference in bivariate latent degradation models. Maker variables are generally chosen
from the available covariates available on a device, making the bivariate model a
type of regression model. The marker is considered to be a ”special” covariate in
that we attribute greater importance to it, because we believe it is more closely
correlated to the degradation variable. At a fixed time point, markers, unlike co-
variates are treated as random variables related to degradation through a parametric
model with unknown parameters. Indexed by time, therefore, a collection of marker
111
variables forms a stochastic process.
Markers, often are taken as functions of several covariates, so called composite
markers. Composite markers can be derived, for example, from principal component
analysis (PCA), factor analysis (FA), or more generally from generalized linear mod-
els. Covariate data scores onto the principal eigenvectors, or latent factors, in PCA
and FA respectively, can be used to represent the composite marker. As another
example, a composite marker can be constructed using the time-dependent cox-
proportional hazards model. In this case, the composite marker variable represents
the instantaneous risk of failure or the hazard rate.
Next we present two less traditional regression structures, that we argue, can
be made available through machine learning methodology.
112
reasonably be approximated by a linear function, like for example degradation or its
drift. It is reasonable, for example, to assume that under higher stress conditions
the degradation drift will increase, but not necessarily linearly.
Consider a cohort of n independent devices tested over a period [0, τ ] as dis-
cussed in chapter 2. The training set D consists of n observations D = [(z i , xi )|i =
1, . . . , n], where z i = (zi1 , zi2 , . . .) denotes the input vector (covariates) and xi =
(xi1 , xi2 , . . .) denotes the output vector for the ith device. In matrix form we can
write D = (Z, x). Note as discussed in chapter 2 the number of input/output ob-
servations made on a device depends on whether it survives or fails in the period
[0, τ ].
x = zT w + (7.1)
where w is the parameter vector of the linear model, f (z) = z T w, and ∼ N (0, σ2 )
is the assumed error distribution between x and f (z) and is assumed i.i.d across
devices. The likelihood of the data given the parameters is given by:
n
Y
L(x|Z, w) = fX (xi |z i , w) = N (Z T w, σ2 ) (7.2)
i=1
113
to the likelihood times the prior, given by:
Zx 1
fw (Z, x) ∼ N , (7.3)
σ−2 ZZ T + Σ−1
z σ−2 ZZ T + Σ−1
z
Z
∗ ∗
fX ∗ (x |z , Z, x) = fX ∗ (x∗ |z ∗ , w)fw (w|Z, x)dw (7.4)
fX ∗ (x∗ |z ∗ , Z, x) ∼ N φT∗ Σz Φ(K + σ2 I)−1 x, φT∗ Σz φ∗ − φT∗ Σz Φ(K + σ2 I)−1 ΦT Σz φ∗
(7.5)
where φ∗ = φ(x∗ ), Φ(Z) is a matrix of columns φ(z). and K = ΦT Σz Φ.
Before we looked at the weight space point of view, which keeps closer with the
linear model perspective. An equivalent way to think about things is the function
space perspective. Here, the weight vector becomes latent and the important concept
is the function itself. So instead of trying to estimate the best posterior w, we
114
Fig. 7.1: Function space view. Left: nine functions drawn at random from a GP prior,
and the dot plots a single observation x. Right: Nine random functions drawn
from the posterior. In both plots the shaded area represents point wise mean
plus and minus two times the standard deviation.
directly estimate the posterior distribution over the functions themselves. This is
possible in GPs where the function is determined by its mean and variance, where
the variance is also known as kernel function and is specified by the user. In the
context of the GP, all random variables, including the r.v. at test data points are
jointly Gaussian, and this means that there are an infinite number of functions
(derived by a specific GP) that can be used to fit the data.
Here we use an example to illustrate the function-space inference process. In
Figure 7.1 observations on x are generated by sampling from a sine function. The
first tier in 7.1 shows the prior and the posterior generated after only one observation
on x, the second tier after 4, and the third after many observations on x. With more
covariate information in the third tier, we observe that the uncertainty in estimation
is reduced.
115
Bayesian linear regression model we can write f (z) = φ(z)T w with prior w ∼
N (0, Σz ). Then
E[f (z)] = φ(z)T E[w]
The covariance function specifies the covariance between pairs of random vari-
ables, for example, cov(f (z p ), f (z q )) = k(z p , z q ) = exp(−1/2|z p − z q |2 ). The spec-
ification of the covariance function implies a distribution over functions. This can
be seen in Figure 7.1, the prior as mentioned earlier is sampled from:
where K is the matrix of kernel functions k(·). The posterior predictive distribution
at test covariates Z ∗ , i.e. the prior conditioned on training covariate observations
Z and new test covariate observations Z ∗ , is given by:
where α = K −1 x.
116
Fig. 7.2: GP mean function m(z ∗ ). The blue dots are the training data (x, z).
117
pendent variable represents the health probability of the device under observation.
The health probability, is an output from the Support Vector Degradation model,
discussed in chapter 10.
Because GPs are computationally very tractable and easy to compute, they
can be used to handle large multivariate covariate observations. Typically covariates
and dependent variables are indexed by time, and therefore, GP regression models
can also be used later for event-time prediction, and we discuss this in more detail
in chapter 10.
Support vector machines (SVM) are based on the idea of large-margin linear
discriminants that seek to find a function f (data) to separate two or more classes
of data by maximizing what is called the hyperplane margin. In this section we
consider linear SVMs, and their possible connection to FHT models.
Consider a cohort of n independent devices tested over a period [0, τ ], and
define training data as: (Z, y) where Z is a collection of r covariate vector ob-
servations on n healthy devices and l covariate vector observations on q failed or
severely degraded devices, such that Z = (z 1 , . . . , z r , z r+1 , . . . , z r+l ). Typically the
healthy training data are collected by observing all n devices early in their life, up
until some predefined time. The degraded training data are collected on all q failed
devices from some predefined time before their observed failure time.
The class label of each z i is given by yi ∈ (+1, −1), where +1 indicates the
membership of z i into the degraded/failed class. The training class label vector y
consists of r healthy covariate observations and l degraded. We assume that Z can
be separated by a decision function f (z; w, b) with appropriate parameters w and
118
b, given by:
m
X
f (z) = (wT z) + b = wi zi + b (7.10)
i=1
subject to
yi (wT z + b) − 1 ≥ 0
n
X
f (z ∗ ) = yi αi z Ti z ∗ + b (7.11)
i=1
where
n
X
αi yi = 0
i=1
and !
n n
1X X
b= yi 1− H ji αj (7.12)
n i=1 j=1
f (z ∗ )
d(z ∗ ) = (7.13)
||w||
119
The drift parameter νX , can therefore, be expressed as a function z ∗ and its
perpendicular distance to the decision boundary f that best linearly separates the
two classes of covariate observations in the training set. One possible configuration
is given by:
νX = β + d(z ∗ ) (7.14)
In equation (7.14) we see that when the distance d(z ∗ ) is close to zero, νX is
not influenced strongly by the covariate observation. For large positive and negative
covariate observations however, νX increases or decreases respectively, and propor-
tionally to d(·). This formulation therefore, incorporates the effects of covariates on
drift, via a non-parametric classification function that discriminates between healthy
and unhealthy/degraded covariate observations. In our opinion, this approach to
including covariates in degradation models is useful in the case of large multivari-
ate covariate data-sets. We anticipate further work on this area as part of future
research.
120
8. VARIABLE THRESHOLD MODEL
Beyond the fixed-threshold model of Whitmore, there is also a need for models
that account for uncertain failure thresholds, which arise in reliability data when fail-
ure is not defined deterministically from the degradation variable. Random thresh-
olds can occur when failure is observed as the result of:
In this chapter, we derive the parametric inference equations for the TMD
data-structure using a variable failure threshold. The model is therefore misspeci-
fied, because it assumes degradation at failure can vary, when in fact its fixed and
equal to a. It is however of practical interest to investigate the efficiency of a variable
threshold model applied to TMD type data-structures. The reason is that failure
thresholds are often arbitrarily chosen, or based on antiquated standards, or more
generally not known. Sometimes engineers can know the range of failure thresholds,
and can prescribe beliefs for the most likely thresholds. In such cases parametric
models for the failure threshold are needed, eg., Gaussian as mentioned earlier.
q p
Y Y
Lθ = fY (S),X(S),T,I(S<τ ) (yi , xi , si , 1) fY (τ ),X(τ ),T,I(S<τ ) (yj , xj , τ, 0) (8.1)
i=1 j=1
122
(y − cx − (νY − cνX )s)2
1
C1 (y, x, s) = p exp −
2πσY2 (1 − ρ2 )s 2σY2 (1 − ρ2 )s
(x − νX s)2
x
C2 (x, s) = p exp − 2
2πσX2 3
s 2σX s
(x − νA )2
1
C3 (x) = p exp −
2πσA2 2σA2
Proof.
fYS ,XS ,T,I(S<τ ) (y, x, s, 1) = fYS ,XS ,S,I(S<τ ) (y, x, s, 1)
fYS ,XS ,S,I(S<τ ) (y, x, s, 1) = fYS ,XS ,S (y, x, s) = fYS |XS ,S (y|x, s)fXS ,S (x, s)
Due to lemmas 4.3 and 4.4, fYS |XS ,S (y|x, s) = fQS (y − cx). Therefore,
For failed items, the level of degradation x at the failure-time s is equal to the
threshold a. Because both X(.) and A are Gaussian, we can replace X(S) with A.
From equation (8.3) we get:
fYS ,XS ,S (y, x, s) = fQS (y − cx)fA,S (x, s) = fQS (y − cx)fS|A (s|x)fA (a)
In the last expression above, fQS (y − cx) is given by term C1 , fS|A by term C2
and fA (a) by C3 .
123
8.2.2 Contribution to likelihood from surviving devices
For each surviving device we observe the pair (y(τ ), x(τ )), and the joint density
for a surviving device is given by:
Z ∞Z τ
fY (τ ),X(τ ),T,I(S<τ ) (y, x, τ, 0) = C4 C5 − C6 C7 C8 dsda (8.4)
−∞ 0
Proof.
fYτ ,Xτ ,T,I(S<τ ) (y, x, τ, 0) = fYτ |Xτ ,I(S<τ ) (y|x, 0)fXτ ,I(S<τ ) (x, 0)
Due to lemma 4.3 the first factor above is given by fQτ (y − cx). Then,
124
X,A A
a0
X(τ)
0 τ time
probability the degradation level belongs to a small interval (x, x + dx) at time τ ,
and that it crossed the failure threshold some time earlier.
125
the threshold variable A.
R∞ Rτ
fXτ ,I(S<τ ) (x, 1) = −∞ 0
fXτ |S,A (x|s, a)fS,A (s, a)dsda
fXτ |S,A (x|s, a) =
= fXτ |S,Xs ,A (x|s, a, a)
= fXτ −Xs |S,Xs ,A (x − a|s, a, a)
X
Due to theorem 4.3 because Xτ − Xs is F≥s measurable it is independent of
X
(S, Xs , A) which is F≤s measurable. Note, as mentioned earlier, the entire process
A(·) is assumed independent of the entire process X(·). Therefore, we can continue
as follows:
Z ∞Z τ
fXτ ,I(S<τ ) (x, 1) = fXτ −Xs (x − a)fS|A (sa)fA (a)dsda (8.7)
−∞ 0
Z ∞Z τ
fXτ ,I(S<τ ) (x, 0) = fXτ (x) − fXτ −Xs (x − a)fS|A (s|a)fA (a)dsda (8.8)
−∞ 0
By plugging in equation (8.8) into equation (8.5) we get the joint density for a
surviving device under a variable failure-threshold. The density fQτ (y − cx) is given
by term C4 , fXτ (x) by C5 , fXτ −Xs (x − a) by C6 , fS|A (s|a) by C7 and fA (a) by C8 in
equation (8.4).
Qq Qp
Lθ = i=1 C1 (yi , xi , si )C2 (xi , si )C3 (xi ) j=1 C4 (yj , xj , sj )×
n R∞ Rτ o
C5 (xj ) − −∞ 0 C6 (xj )C7 (xj )C8 dsda
126
estimates of the process parameters θ̂.
127
9. SUPPORT VECTOR DEGRADATION MODEL
9.1 Introduction
129
The Bayesian SVM algorithm is trained in the absence of failure data (negative class
data), as is the case in many mission-critical systems. The contribution of this work
to the field of reliability are the following:
130
in X. The dimension for each model can be chosen to be as low as 1, and as high
as m − 1, with m potential models for each, respectively, considering all the possible
combinations. Here, for expositional simplicity, only two models are considered, [M ]
and [R], each of which is two-dimensional.
PCA was chosen to preprocess the original data to extract features related
to changes in the variance of the data. In the context of statistical control theory,
the variance and its changes are strong features indicating the onset of anomalies in
multivariate systems [82]. Other options for this step include blind source separation
(BSS), and independent component analysis (ICA), or, more generally, generalized
linear models (GLMs). PCA is a special case of GLM, and although it suffers from
the assumption of linearity and normality of the data (situations that are arguably
not often encountered in real data sets), transformations can apply to the original
data to approximate normality [83].
When failure data (negative class) are not available, a kernel density estimate
(KDE) is computed for the projected (positive class) training data in the two sub-
spaces [M ], and [R] to estimate the likelihood of the projected data, and from it
to construct the negative class. The SV classifier constructs two predictor mod-
els, D1 and D2 , for each subspace. A soft decision boundary is constructed by
fitting the training data with a model for posterior class probabilities using a lo-
gistic distribution that maps classified data to posterior classification probabilities:
PM , PR1 , . . . , PRM , respectively. The joint class probability Jp from the subspaces is
used for the decision classification.
Support vectors produce an uncalibrated value that is not a probability. There-
fore, the algorithm uses the support vector decision function, D, to produce a
posterior probability, P (class|input), according to a Bayesian formulation. Fi-
nally, the joint posterior class probability can be weighted with a weight vector
W = [w1 , . . . , wm ] to emphasize some models as opposed to others. This weighting
131
M pM
w1
R1 pR1 w2
w3
X R2 pR2 Jp
wm
Rm pRM
Fig. 9.1: Algorithm flow diagram showing the processing of the data.
could be beneficial for emphasizing the results from models, usually the principal
model [M ], which captures more of the data covariance information. In this paper,
all models are weighted equally (W = I).
In this chapter, we decompose the training data into two lower dimensional
subspace models: the principal model [M ], and the residual model [R]. We use
singular value decomposition (SVD) of the input data, X [84], [85], [86], [87]. The
SVD of data matrix X is expressed as X = U SV T , where S = diag(s1 , . . . , sd ) ∈
Rn×d , and s1 > s2 > . . . > sd are the ordered singular values. The two orthogonal
matrices U , and V are called the left, and right eigen matrices of X. Based on the
SVD, the subspace decomposition of X is expressed as:
The diagonal matrix, SM , are the singular values (s1 , . . . , sk ), and (sk+1 , . . . , sd )
belonging to the diagonals of SR . Any vector X can then be represented by a
summation of two projection vectors as shown in equation (9.2), where PM = U U T ,
and PR = I − U U T are the projection matrices for the principal, and residual model
132
subspaces, respectively. Both subspaces comprise the total data dimension. In this
framework, we can apply a SV classifier, having oriented the data such that we
can better capture system failures that are reflected in changes in variance, and so
that we can ”break down” the effects of multivariate data into separate models,
each examining a different effect of the data on changes in variance. Then we can
envision combining the results in the end to achieve a ”global” detection result, as
we demonstrate later in this paper.
X = PM × X + (I − PR ) × X (9.2)
SVMs alleviate the need for algorithms with statistically grounded frameworks,
algorithms that require knowledge of the distribution of the random variables. SVMs
are based on the idea of large-margin linear discriminants that seek optimum margin
hyperplanes where the separating plane is chosen to minimize a risk bound moti-
vated by structural risk minimization. Nonlinear extensions were introduced by the
authors in [88], and [89] with a generalization often referred to as the ”kernel trick”,
which builds on a direct consequence of Hilbert’s space theory. Here, we review the
linear SVMs to highlight certain concepts that we will use in this chapter. Given
the data structure defined earlier, linear SVMs apply a linear model f that maps a
m-dimensional real valued vector to a binary scalar as shown in equation (9.3).
T
f (x) = sign(wSV M x + (x)) (9.3)
n
X
wSV M = αi yi xi (9.4)
i=0
T
In (9.3), wSV M x = D(x), and its weights wSV M are given by equation (9.4),
are normal to D; b/||w|| is the perpendicular distance from D to the origin, and
133
||wSV M || is the Euclidean norm of wSV M . The margin in classification is the distance
between the nearest positive and negative labeled data points. For linearly separable
cases, training the SVM is performed by solving the following optimization problem:
argmin ||wSV M ||, and the constraints are combined into a set of inequalities ∀ i:
yi (xi w + b) − 1 ≥ 0.
Training the SVM becomes an optimization problem given by equation (9.5),
and constrained by equation (9.6).
n
X
2
(α1 , . . . , αn , b) = argmin(1/2||w|| + C ξi ) (9.5)
i=1
yi (wT xi + b) + ξi ≥ 1 (9.6)
Lastly, in the case where a linear decision function is not suitable for the data,
the above methods can be generalized using a transformation to another Euclidean
space using a map function called Φ, where the training data are linearly separable.
More reviews of SVMs can be found in [90], [91], [92], and [93].
Evidence Framework
134
normally distributed around zero with a certain variance, therefore modeling the
conditional density of Y |X also as normal. c) Because of assumption b), we can
express the posterior class density as a function of an SVM related term, namely ξi .
To see this result, we consider the training data Z, assume that the joint P (Z, w)
exists, and assume that the conditional P (w|Z) can be expressed as
P (Z|w)P (w)
P (w|Z) = ∝ P (Z|w)P (w) =
P (Z)
= P (y|X, w)P (X, w)P (w) (9.7)
The posterior on the weights can be expressed as the product of three dis-
tributions, as shown above. The probability density over observations given the
parameters is modeled through a binomial distribution to account for the possible
states of the response random variable Y , and is given by (9.8) with 0≤ q(x, w) ≤ q1,
and q(x, w) = P rob(y = +1|X, w).
n 1 + yi 1 − yi
Y
P (y|X, w) = q(xi , w) 2 (1 − q(xi ), w)) 2 (9.8)
i=1
135
tive function for an SVM, and is given by (9.11).
n φ 2
Y ||w|| )
(−
P (w|Z) ∝ (ρa (1 − ρ)b )e 2 (9.10)
i=1
n
X n
X
−log(P (w|Z)) =C − alogρ − blog(1 − ρ) (9.11)
i=1 i=1
ρ = erf c(−wT xi )
φ
C= ||w||2
2
1 + yi
a=
2
1 − yi
b=
2
By considering the asymptotic expansions for the error functions above, and if
it can be shown that the expansions of the two sums reduce to a function of ξi , the
log posterior on the weights of the linear model has an equivalent form to that of
the SVM optimization in equation (9.5). This connection is useful because it effec-
tively lays down a strong informative prior for modeling posterior class probabilities
of future test data. This prior is implemented by treating D(x) as the optimum
classifier for the given training data. This fact will be used later in the chapter to
provide rationale for the design of a posterior classification probability given D(x).
In many real world systems, especially mission critical systems, and compo-
nents for which failures are not known, training data consists only of the positive
class. To obtain estimates for the failure space (negative class), novelty detection
136
as discussed in [98], [99], [100], [101], and [102] among others (see [103], and [104]
for more general review) is approached primarily as a data-versus-density problem,
where the negative data are assumed to be generated from an unknown distribu-
tion, say Q(X). The density of Q is intentionally left uninformative, and usually
uniformly distributed to reflect the lack of any prior knowledge about anomalies.
Authors in references such as [98], and [102] discuss sampling schemes for Q
that optimize supervised function estimation techniques (e.g., SVM) to best infer
a general classification boundary for the given positive class training data. The
sampling approaches depend on the choice of a prior for Q, which in the absence of
any evidence is measured on the entire metric space spanned by the positive class
training data, and suffers from high dimensionality. In [102], the authors discuss a
negative class selection algorithm for data collected from various Internet sites. In
this work, unlabeled data was made available by sampling the Internet, which is
different from the situation we describe here. In this paper, there are no unlabeled
data, and we cannot sample from a universal set (the Internet, for example).
Other approaches, as mentioned earlier, use the origin as the negative class in
the applied feature space induced by some kernel function [105]. Others [106] extend
this idea, and assume that all data points close enough to the origin are also con-
sidered as candidates for the negative class. Some of the critiques, however, of the
one-class classification approach motivated by [105] focus on its sensitivity to specific
choices of representation and kernel in ways that are not very transparent [106]. Fur-
ther, its assumed homogeneous input feature space relies on comparable distances
between data, which can lead to inaccurate classifications with non-Gaussian dis-
tributed data [83]. The authors of [83] propose a rescaling of the data in the kernel
feature space to make it robust against large-scale differences in scaling of the input
data. The data are rescaled such that the variances of the data are equal in all
directions using kernel PCA.
137
The primary approach to one-class classification has been largely based on the
work discussed above. In essence, the problem reduces to making the most of the
information at hand, the positive class training data. As such, it becomes impor-
tant to extract features from these data that can improve inference about potential
anomalies. An important feature of the training data is its density, which can be
estimated computationally, although this is an expensive task in high dimensions.
Therefore, on a practical level, the estimate of the negative class is seen as a
conservative representation of a potential system failure space, an assumption that
could lead to poor generalization of the algorithm in situations where the predictor
model is not updated to reflect changes in the system performance characteristics.
Such changes are plausible, for example, in a reliability setting in which the system
has aged so that its performance signature has changed, but it is still functioning in a
”healthy” state. Another example is a case where the original training data were not
complete enough to represent the global system performance regimes (universal set),
and in such situations the predictor model will naturally fall victim to large numbers
of false alarms. Therefore, a one-class-classifier approach to novelty detection must
be subject to complete, updated training data.
To utilize SVMs for classification, the negative class must be estimated first
by considering the density of the positive class (training data) following similar
reasoning as the authors of [107]. This work can be accomplished in several ways,
one of which is to use a kernel density estimate (KDE) of the training data through
the use of Gaussian kernel functions. For this work, the negative class was estimated
based on assumptions on the failure space, summarized in 9.1
Definition 9.1. The failure space is a) not linearly separable from the healthy train-
ing data, b) prevalent in the space not occupied by the healthy training data, and
therefore c) assumed to conform to the distribution of the healthy training data.
Through this definition, we aim to achieve minimum volume sets (mvs), similar
138
to work discussed in [108], [109], [110], and [105], that find sets of density functions
that correspond to regions with the minimum volume or Lebesgue measure for a
given error 1 − α [109]. An mvs in a class of measurable sets C for an error α is
defined in reference [110] as
139
influence weighted according to their Euclidean distance from x. One choice for a
smooth φ is the standard normal distribution N (0, 1).
To overcome over-parameterized density estimates that do not generalize well,
the bandwidth, h, is determined through a nearest neighbor approach in which h is
√
selected as the value that produces a volume around xi containing n neighbors.
This approach personalizes the value of h to each data point xi , and effectively
smoothes out the density in areas with sparse training data information.
The negative class data are constructed by selecting grid coordinates where
the likelihood ratio ρ of the training data is below a threshold, τ , which is a grid
center, and is labeled as a member of the negative class if ρ ≥ τ . The likelihood
ratio is the ratio of negative to positive posterior class probabilities, as shown in
(9.13). The denominator is computed by the KDE, and the numerator is modeled
as a function of the gradient of the likelihood function.
In (9.13), we use P (Y = −1|X = xi ) as p−
i , and P (Y = +1|X = xi ) as
p+
i . We note that the model favors the numerator proportional to the square of the
likelihood function gradient, ∇L. Practically, this means that, in areas where the
likelihood function changes faster, the negative class is more similar to the positive
class.
P (Y = −1|X = x)
ρ= (9.13)
P (Y = +1|X = x)
p− + + 2
i = (1 − pi ) + (1 − pi )∇ L (9.14)
Once the D(x) is constructed through an SVM using the positive and esti-
mated negative training data, the argument is, as motivated earlier by its statistical
140
properties, that D(x) is the optimal classifier.
−1 if P (Y = +1|X = x) < 0.5
=
+1 if P (Y = +1|X = x) ≥ 0.5
141
D(x)
Fig. 9.2: Logistic distribution model for posterior class probabilities.
total probability, re-expressing the sum in the denominator, we get a function with
parameter β that is a logistic-type distribution.
P (X = xi |Y = +1)P (Y = +1)
P (Y = +1|X = xi ) = P
P (X = xi |Y = a)P (Y = a)
a=−1,+1 (9.16)
1
=
(1 + exp(βi ))
P (X = xi |Y = +1)P (Y = +1)
βi = log (9.17)
P (X = xi |Y = −1)P (Y = −1)
142
1
pi = P (Y = +1|X = xi ) =
1 + e(−a1 g(xi )+a2 )
n
Y
P (c1 , . . . , ck ) = pci i (1 − pi )1−ci (9.19)
i=1
143
According to Bayes’ rule, the conditional class probability is given by
P (XM , XR |Y )P (Y )
P (Y |XM , XR ) = =
P (XM , XR )
P (XM |Y )P (XR |Y )P (Y )
= P
P (XM , XR |Y )P (Y )
y={−1,+1}
P (Y |XM )P (Y |XR )P (Y )
=P (9.21)
P (Y |XM )P (Y |XR )P (Y )
y
144
dation. In section 7.3 and 7.3.1 we consider SVMs as a potential regression structure
in FHT models, here we argue that the posterior class probability (a direct result
of SVM classification) can be used as a marker.
145
10. CASE STUDIES
0.80.8
0.60.6
Probability
0.40.4
0.20.2
0.00.0 0 25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 500 525
800 900 1000 1100 1200 1300
800 900 1000 1100 1200Observation
1300
#
Observation #
Fig. 10.1: Joint posterior class probability vs. observation for Lockheed Martin test data
set
CALCEsvm LibSVM
1
100.0%
0.8
80.0%
Probability
0.6
Probability
60.0%
0.4
40.0%
0.3
20.0%
0.0
0.0%
0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 320
0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 320
Observation #
Fig. 10.2: Joint posterior class probabilities for CALCEsvm, and the open source support
vector classification software called LibSVM.
CALCEsvm was used on these data. Fig. 10.1 shows the detection results.
The algorithm detected the first two periods of anomalies, namely those between
912 and 1040, and between 1092 and 1106.
CALCEsvm was compared to the open source support vector classification
software called LibSVM [93]. The setup for LibSVM used its two-class C-SVC
147
setting with input, the training data used in CALCEsvm. Because the one-class
SVM in LibSVM does not provide posterior class probabilities for test data, we
compared the two-class classification between CALCEsvm and LibSVM. For this
comparison, the negative class training data were taken from the output estimate
of CALCEsvm, and used as the negative class in LibSVM. Therefore, the actual
comparison was made between the two-class SVM algorithms of CALCEsvm, and
LibSVM. The option settings used for LibSVM are listed below with the margin
penalty parameter, and tolerance setting of the termination criterion parameter
chosen arbitrarily, and kept the same for both CALCEsvm and LibSVM.
The accuracy comparison was performed through three tests: 1) a direct com-
parison of the quadratic optimization results: the objective function, the sum of
the Lagrange multipliers, and the number of support vectors; 2) detection accu-
racy based on class index only; and 3) detection accuracy based on the range of
probabilities.
In Table 10.1, b0 is the bias term, w2 is the objective function equal to αT Hα
where α ∈ R1×n is the Lagrange multiplier vector, and H ∈ Rn×n is the Hessian
matrix, where n is the length of the SVM training data. The parameter is the
tolerance of the termination criterion, and nSV is the total number of support
vectors. The results in table 10.1 show that the performance of the software is
148
Tab. 10.1: SVM Optimization Results
Principal Model Residual Model
CALCEsvm LibSVM CALCEsvm LibSVM
b0 0.181 181 0.099 0.099
w2 15.40 7.70 14.70 5.60
0.0001 0.0006 0.0001 0.0005
nSV 24 23 14 14
comparable to the difference found in the objective function. The number of support
vectors, and the bias term were found to be the same.
The second, and third tests compared their detection accuracy against the
known periods of anomaly. Each test file was coded with a column variable z ∈
{−1, +1}, indicating the known class of each observation, an index of +1 for the
healthy data, and -1 for the anomalous data. LibSVM counted the number of
misclassified observations based on the coded variable z. Table 10.2 shows the results
comparing LibSVM to the CALCEsvm output. The first column in the table shows
the detection accuracy based only on the class index, whereas the second column
shows the detection accuracy based on a probability index. In the first comparison,
both performed almost identically (see first column in table 10.2), but the second
comparison (second column) clearly favors CALCEsvm. This result can be seen by
comparing the accuracy of 98.1% for CALCEsvm vs. 30.5% for LibSVM given the
criteria that the posterior class probability for a test observation should lie within
the range specified in the algorithm, here 0.8 to 1.
The second comparison was performed based on a probability index reflecting
an ”expert” knowledge of system ”health”. This index therefore pertains to a belief,
and is subjective to the user. Nonetheless, this index is based on an intuitive argu-
ment: because the posterior class probabilities reflect the certainty/uncertainty of
the classification/detection, a known ”healthy”, and or known ”unhealthy” observa-
tion should be associated with high, and low probabilities (or ranges of probabilities),
149
Tab. 10.2: Comparison of CALCEsvm and LibSVM Detection Accuracy Against Lockheed
Data
Detection Accuracy Results Based On
Class Index Probability Probability Index Range
CALCEsvm
100.0% 98.1% 0.8-1.0
100.0% 0.0-0.4
LibSVM
99.6% 30.5% 0.8-1.0
100.0% 0.0-0.4
respectively.
In the Lockheed data set, there are two system levels: ”healthy” and ”failed”.
Both system levels are known for the whole data set. The ”healthy” level is set to
be represented by posterior class probabilities between 0.8 and 1, and the anoma-
lous level by probabilities between 0 and 0.4. Stronger restrictions can be modeled
by expanding the range for the anomalous level, and shrinking the range for the
”healthy” level. In light of these explanations, CALCEsvm had a 1.9% error rate in
its detection accuracy as opposed to 69.5% for LibSVM. The reason LibSVM per-
formed at 30.5% accuracy is because two out of three periods with ”healthy” level
operation were captured (by LibSVM) with a posterior class probability at around
0.75 to 0.78, therefore falling short of the user-defined ”healthy” range of 0.8 to 1.0,
and failing to correctly classify the healthy periods. LibSVM, as did CALCEsvm,
captured the failed periods with 100% accuracy.
A second case study was performed using simulated correlated data consist-
ing of three random variables from three different but s-dependent distributions to
construct the training data set. The objective in this case study was to test the
algorithms on a system that was degrading, and in which the degradation took
150
place in the presence of considerable noise. Copulas were used to build a simulation
model consisting of three random variables: Gamma(2, 1), Beta(2, 2), and t(5). The
family of bivariate Gaussian copulas is parameterized by ρ = [1 ρ; ρ1], the linear
correlation matrix. The random variables U1 , and U2 approach linear s-dependence
as ρ approaches +/-1, and approach complete statistically independence as ρ ap-
proaches zero. The Gaussian, and t copulas are known as elliptical copulas, and
can generalize higher numbers of dimensions. Here we simulate data from a trivari-
ate distribution with Gamma(2, 1), Beta(2, 2), and t(5) marginals using a Gaussian
copula.
Test data were generated from the trivariate distribution of Gamma, Beta,
and t random variables; and were set up such that three degradation periods were
generated. The first period was designed to be ”healthy”, the second introduced
a shift in the mean for each variable separately while maintaining the correlation
structure, and the third period introduced a larger shift in the mean.
The CALCEsvm results are shown in Fig. 10.3, with the four periods identified
by breaking perforated lines and an index P 1 through P 4, where P 1 is the identifier
for the ”healthy” period, with mean equal to nominal, and P 2 through P 4 having
successively increasing changes in the mean.
The results of the algorithm show the ability to capture the trend of simulated
degradation in the presence of noise. The beginning period that shows a dip in the
probability estimate is a direct result of an initial over-smoothing (implementation
of the exponential smoothing), and can be ignored for practical purposes. The larger
result is the algorithm’s ability to correctly classify the data for each period of oper-
ation, and to capture the expected trend. CALCEsvm results were compared to the
results obtained from LibSVM, and are tabulated in table 10.3. The probabilities,
as in the Lockheed Martin case study, again reflect a belief about the interpretation
of the posterior class probabilities. In this case, posterior class probabilities between
151
1
0.9
0.8
0.7
Probability
0.6
0.5
0.4
0.3 P1 P2 P3 P4
0.2
0.1
0
0 50 100 150 200 250 300 350 400
Observation number
Fig. 10.3: Joint positive posterior class probability for simulated data set.
0.8 and 1 are acceptable for a ”healthy” system, probabilities between 0.7 and 0.85
are acceptable for the next level of ”health” allowing for some overlap, and so on
until the range between 0 and, say 0.5 for example, are used to classify the system
as failed.
The comparison of accuracy results based only on the class index shows that
both algorithms performed virtually identically for the given probability ranges.
Both CALCEsvm, and LibSVM had a detection accuracy rate of 100% in P 1; in
P 2, both algorithms performed noticeably poorly; and both improved in P 3 and
P 4 to 88% when the degradation became more distinct. A comparison of accuracy
results based on the posterior class probabilities shows a slight improvement in the
152
Tab. 10.4: CALCEsvm Accuracy Results for Simulated Data
CALCEsvm Output Based on
Class Index Probability Probability Index Range
100.0% 81.2% 0.8 - 1.0
3.8% 31.7% 0.7 - 0.85
29.0% 31.0% 0.3 - 0.7
88.0% 84.0% 0.0 - 0.4
LibSVM CALCESVM
98.0%
94.0%
Accuracy
90.0%
86.0%
82.0%
0.5 0.55 0.60 0.65 0.70 0.75 0.8
Start Range
Fig. 10.4: Joint positive posterior class probability for simulated data set in P2.
performance of each algorithm for P 1, and about the same performance for the
other periods. Fig. 10.4 plots the detection accuracy of CALCEsvm and LibSVM
vs. the start value for the probability index for levels 1, 2, and 3. For example,
from the plot, it can be seen that when the lower bound on the probability index
is 55%, and the upper limit is fixed at 100%, CALCEsvm has a detection accuracy
of 96% vs. approximately 89% for LibSVM. This is a very liberal bound, as it says
that any posterior class probability above 55% can be used to classify a test point
as ”healthy” instead of anomalous to some degree.
In Fig. 10.4, the x-axis shows the varying lower bound for period 2 (P 2), and in
Fig. 10.5 the varying lower bound for period 3 (P 3). The y-axis shows the accuracy
of the algorithms in classifying test data. As the lower bound on our belief becomes
153
100%
80%
Accuracy
60%
40%
20%
0%
0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80
Start range
Fig. 10.5: Joint positive posterior class probability for simulated data set in P3.
more stringent (that is, we require higher certainty in the prediction), the accuracy
of the algorithms falls. Because the anomalies are more distinct (due to stronger
outliers, and reflected by lower posterior class probability values) in P 2 than in P 3,
as the lower bound is ”tightened” similarly for both periods, both CALCEsvm, and
LibSVM perform better in P 3. Here, for example, when the lower bound on the
probability index was 70%, and the upper held at 100%, CALCEsvm had lower than
94% detection accuracy, whereas LibSVM had an accuracy of 89%.
The last set of data that we used to test and validate the BSVM algorithm
came from the NASA data repository of degradation data. The data was and is still
available as part of a competition. The data is composed of training data and test
data. The training data consists of multivariate time series of covariate observations
from different enginesa total of 218 engines, which we call units. Each engine de-
graded due to wear based on the usage pattern of the engines and not necessarily due
154
to any particular fault mode. Each unit started with different unknown degrees of
initial wear and failed at an unknown level of wear. Noise was injected into the data
to represent manufacturing variation, process noise, and measurement noise. The
engines were exposed to six operational settings that provided information about
the operational mode and environmental conditions of the engine. The objective
was to predict remaining life given a sequence of covariate measurements for a test
engine.
Modeling degradation using a classification approach (as discussed here) we
need to address some potential limitations (lessons learned). The classifier becomes
inaccurate: 1) if the variance of the initial wear is large, and 2) if the difference in
performance measurements between healthy and unhealthy engines is not significant.
This is especially true with this algorithm, which is based on a binary classifier, i.e.,
one that discriminates on the basis of two classes. The data was first grouped
according to the operational setting number 3 into six sets, each representing data
collected at each setting respectively. Setting number 3 takes six possible levels: 0,
20, 40, 60, 80, 100 and represents some type of stress condition. To account for
the six settings (stress levels), six classifiers were used based on the training data
obtained from each group respectively, and each units level of health was estimated
for each setting separately. For each setting the healthy training set was taken as a
percentage of the initial observations, and the unhealthy/degraded data were taken
as a percentage from the final observations across all units. The test data were
similarly partitioned.
To get predictions, we first trained a classifier based on CALCEsvm. We used
the output of CALCEsvm, the joint posterior class probability, to make predictions.
As discussed above, CALCEsvm reduces the dimensionality and outputs a univariate
time series of health estimates. In this case, CALCEsvm was used as a two-class
classifier, using both healthy and unhealthy training data. Figure 10.6 plots the
155
Fig. 10.6: Distribution of projected multivariate data on to first two principal components
distribution of the healthy, unhealthy, and test data for setting 0 projected onto the
first two principal components. In Figure 10.6 we see that the two sets are suitable
for classification-based detection approaches. We also see the trajectory of the test
data, starting in the healthy set and moving towards the unhealthy. This pattern
presents a time series of health estimates that exhibits a trend which is suitable for
prediction.
CALCEsvm was applied to the training data for each setting and gave a se-
quence of health estimates indexed by a unique cycle number associated with each
observation. By combining the health estimates from each setting into one vector
and sorting by the cycle number, we got the final time series of health estimates
for each unit. For example, figure 10.7 plots the probability of health for each test
observation across all settings for training unit 200. In figure 10.7 we get the desired
and expected drift in the health estimate as the unit ages. Health estimates close
to 0 indicate a failing or failed unit.
A similar time series was estimated for each test unit from 1 to 218. There
156
Health Probability
Cycle #
were some test units with very short histories that had not started to degrade in
health by last observation. For these units we did not get clear downward trends,
and prediction in these cases was difficult. To predict the expected time to failure,
we fit a GP model to the resulting time series. To account for health estimate
variability, we fit the GP multiple times, each time perturbing the estimates with
Gaussian noise with a mean and variance estimated from the training data across
all units (cross unit variation). The cross unit variation for all test units is shown
in figure 10.8. We see that the variance increases towards the end of life, validating
the simulation design [38]. The drop in the expected variance, as seen in the lower
plot of figure 10.8, is a result of a decreasing sample size during those time points
(some test units survived longer than others).
Similar results were obtained for the training units. Using this variation profile,
we injected noise into the original health estimates and fit the GP model multiple
times. Each fit will give us the mean and the 2 standard deviation paths. We mod-
eled the failure-time as the first time the mean path hit some fixed threshold; for
example, a health level of 0. Figure 10.10 plots an example of such a procedure; it
shows the GP fit to the health estimates for unit 200. It also overlays the empirical
distribution of the 50 first hitting times to a threshold level set at 0. The expected
157
1
0.8
Probability of Health
0.6
0.4
0.2
0
0 50 100 150 200 250 300 350 400
Cycles
Fig. 10.8: Variation in estimated health probability across all training units
0.35
0.3
0.25
Variance
0.2
0.15
0.1
0.05
0
0 50 100 150 200 250 300 350 400
Cycles
158
0.9 0.05
0.045
0.8
0.04
0.7 0.035
0.03
0.6
Probability of Health
0.025
0.5
0.02
0.4 0.015
0.01
0.3
0.005
0.2 0
250 300 350 400
0.1
-0.1
0 100 200 300 400 500 600 700 800 900
Cycles
Fig. 10.10: GP model fit to the health estimate time series of a test unit
159
11. MULTISTATE MODELS AS DEGRADATION MODELS
11.1 Introduction
161
probability of an device given its marker process up until time t. Doksum and Hoy-
land (1992) [116], Doksum and Normand (1995) [117], Whitmore (1995) [118], and
Whitmore et al (1998) [24], among others, model degradation by a Wiener diffusion
process. Lee, Degruttola and Schoenfel (2000) [26], Henderson (2000) [120] consider
extensions to the bivariate Wiener diffusion models with time-varying longitudinal
marker processes. Longitudinal data are collected for each marker process not only
at failure but also throughout the lifetime of the device. Satten and Longini (1996)
and Hendriks (1996) use Markov models to combine a longitudinal trajectory with a
survival event. In this work the authors develop parametric and predictive inference
models that are generally analytically solvable. They make predictions of time to
failure, expected time to failure, and other functionals of time by using the maxi-
mum likelihood estimates (MLEs) of the model parameters and plugging them into
the predictive equations.
Meeker and Escobar (1998) [121], Commenges (1999) [122], Commenges (2002)
[123], Bagdonavicius et al (2002) [124], Putter et al. (2006) [125], Pena (2006) [126],
Machado et al. (2008) [128], Andersen and Perme (2008) [129], Cook, Lawless,
Lakhal-Chaieb and Lee (2009) [130], Aalen, Borgan and Gjessing (2008) [131], Cook
and Lawless (2007) [132] consider multistate models for survival and event history
analysis based on counting processes. An excellent exposition, review and applica-
tion of counting processes are given by Andersen, Borgan, Gill and Keiding (1993)
[133], and Aalen and Johansen (1978) [134]. The counting process framework pro-
vides a non-parametric approach to inference and prediction, most famously through
the Kaplan Meier (KM), Nelson Aalen and Aalen-Johansen (AJ) estimators. The
AJ estimator or product-integral as it is otherwise known, estimates the transition
probability matrix of a nonhomogeneous Markov process. This theory is useful to us
in developing a simple markov multistate model to compute 1) the expected time-
to-failure, and 2) the survival probability, for a surviving device given its marker
162
process up until time t.
The multistate model is a natural and simplifying representation of the bivari-
ate marker-degradation process. With the multistate model, the multivariate state
space of the marker and degradation processes is re-parameterized to a univariate
Markov chain by expressing the state space as a list of all possible points of the
marker-degradation vector pairs. Longitudinal marker observations are naturally
incorporated into the model by directly contributing to the estimate the transition
probability matrix, making predictions dependent on the process history. The time
to failure is modeled as the FHT of the degradation variable to the failure threshold.
As discussed earlier, because degradation variables that define failure are typ-
ically not observable in computer systems, we examine the situation in which the
degradation variable is unknown and unobservable. This might be the case in most
complex electronic systems, such as computers, that can exhibit intermittent events
indicating failure, but the underlying mechanism and therefore degradation variable
is unknown. In this case predictive inference is based only on marker observations.
In this work inference and prediction are based on data from one computer, and each
failure time is modeled as independent and identically distributed. In our proposed
model we assume that the system is as good as new after each failure. Critical error
messages generated by internal performance monitoring software are considered as
intermittent failure events correlated to failure.
Traditionally failure time models do not accommodate dynamic model pa-
rameters. In other words, most failure time models assume unknown but fixed
parameters that are not influenced by changing environments. The focus in current
literature is to incorporate information about the environment and or about the
performance level of the system or device under observation, in order to model time
dependent model parameters. Typically separate models are embedded into the
time to failure probability density function to express the dependence on covariates
163
(sensor information) or markers like is the case in this study. Validating models and
assumptions in such approaches becomes important and can be a limitation.
Computer failures are defined from a user perspective, in which failure is seen
as the loss of functionality of the computer. Failures are defined as automatic or
forced restarts of the system. Automatic restart is triggered by the computer while
forced restart is performed by the user. The failure time is defined as the time at
which one of these events occurs. Once the computer is restarted, the process starts
from time zero again, until the next failure. The time to the next failure can also
be considered the time between failures.
The hypothesis is that each computer will fail (as defined above) as a result
of usage and environment stresses. The hypothesis, therefore, is that there exists a
measurable variable whose values correlate to the failure time, in other words, that
there exists a marker/precursor to failure. Error messages generated by internal
monitoring software are believed to be precursors to failure in computers. Fault
events are simply called errors in the remainder of the chapter.
Lifetime and covariate data are only collected from one computer system. As
mentioned earlier, the computer is assumed to be as good as new after each failure,
and therefore we also assume that each failure-time is independent and identically
distributed. For expositional simplicity, in this chapter, we use the term ”device” to
represent the information associated with the computer between each failure. For
example, device 1 represents the computer between time zero and the time of the
first failure. device 2 represents the computer from right after the first failure up to
the second failure, and so on. Note also that this formulation is valid in the case
when we have failure-time and covariate data from a sample of computer systems,
in which case the assumption of independence is stronger, and the modeling holds
164
the same.
165
Y
X
a
Fig. 11.1: Illustration of sample paths from a bivariate stochastic process {X(t), Y (t)}
rence of a future problem. For example, a Warning event message is logged when
disk space starts to run low. c) Error, an event that describes a significant problem,
such as the failure of a critical task. Error events may involve data loss or loss of
functionality. For example, an Error event is logged if a service fails to load during
startup. Examples of error events are: Application Hang, Memory Access Denied,
etc.
A FHT model can describe the relationship between marker, degradation and
lifetime variables. Typically, lifetimes are modeled by the FHT of the degradation
process to a given threshold level a, namely T = min(t|X(t) = a), in other words the
minimum time at which the degradation process reaches the threshold level. Typi-
cally, the fht model for a joint stochastic process {X(t), Y (t)} (degradation/marker)
needs to condition on observations on Y (t). Each stochastic process is defined on a
state space S and time space T . The state space of the marker variable is the set
of all natural numbers, which record the cumulative number of errors.
Figure 11.1 illustrates the bivariate stochastic process X(t), Y (t) with a
166
p00(t) p11(t) p22(t) p33(t)
p01(t) p12(t) p23(t) p34(t)
0 1 2 3 ... k
p0T(t) pkT(t)
T
threshold level for X(t) at a, and two possible paths, one resulting in a failed device
and the other a surviving device. Because the state space of Y (t) takes discrete
states, the bivariate process looks like a step function. Later, we present a case
study where we generate observations on Y (t) by modeling it as a Poisson process.
The discrete state space of Y (t) can be expressed as a multistate representation as
illustrated in figure 11.2.
Each state in the multistate model represents the cumulative number of er-
rors experienced by the device, and the devices latent degradation level, with the
exception of the last state which represents the cumulative number of errors and
a degradation level of a. At any moment in time, therefore, the multistate model
accounts for two possible events; the occurrence of an error or a failure. There are
three possible transitions at any moment in time: a) to the failure state F , b) to the
right adjacent state from 0 to 1 or from 1 to 2, etc., and c) remain in the same state,
i.e., via p11 (t) from state ”1” to state ”1”. The characteristic probability transition
matrix for this model is a sparse banded symmetric matrix P with elements along
the primary and upper adjacent diagonal and along the last column of the matrix
as illustrated in figure 11.3. Further description of the resulting matrix is given in
section 12.6.
Although in principle there may not be a maximum number of states for the
167
marker variable, in this model a maximum number of states is used based on the
sample data. For example, k in figure 11.2 is equal to the maximum number of
errors across all devices for the computer under test. We assume that all devices
begin in state 0 at t = 0, that state 0, . . . , k are transient, and that state F is
absorbing. With the multistate model, the bivariate state space {X(t), Y (t)} is
reparametrized to a univariate Markov chain by expressing the state space as a list
of all possible points in the bivariate space. In other words, each state represents
the cumulative number of observed errors and the associated state of health. A
univariate representation of the state space allows for a model to treat the latent
state of health (degradation) abstractly, without making any assumptions about the
joint distribution of degradation and marker, and make predictions therefore about
degradation purely through the marker.
Consider a Markov chain Yn for the multistate model, with state space Ω =
{ω ∈ N } and transition probability matrix P = {pij } i, j ∈ Ω, with pij ≤ 0, Σk∈Ω
pik = 1 for all states i, j, pij = P (Yn+1 = j|Yn = i). The transition probability ma-
trix is a sparse banded matrix with elements along the primary and upper adjacent
diagonals and along the last column of the matrix as illustrated in Fig. 11.3.
Each element in the transition probability matrix represents the probability
of transitioning into the corresponding state in the multistate model. Elements
in row i represent states j 6= i that the system transitions to from state i. For
example, if we consider row i = 1, then element (1,2) of the matrix represents the
probability of transitioning from state 1 to state 2, and element (1,1) the probability
of transitioning from state 1 to state 1, in other words remaining in the same state.
The matrix is sparse because most of it is populated with zero elements. Zero entries
in the transition matrix are used to model a zero probability of transition into the
168
corresponding state. For example, in our model, element (1,3) is a zero element
because the system cannot transition from state 1 to state 3; physically it cannot
experience the third error before it experiences the second error. The diagonal
entries then are the probabilities of transitioning along the marker states and the
last column the probabilities of failure from each respective state. The multistate
model is also an absorbing Markov chain.
where E is a matrix written in terms of Q and U , but its expression is not needed
here. The form P n shows that the entries of Qn give the probabilities of being in
each transient state after n steps.
169
p11 p12 0 p1n
p22 p23
0
0 0 p33 .
. .
. .
.
pn-1n
0 0 pnn
Theorem 11.1. In an absorbing Markov chain, the probability that the process will
be absorbed is 1 (i.e., Qn → 0 as n → ∞).
The results from equation (11.2) concern the probability of remaining forever
in the transient set, or alternatively, the probability of never being absorbed by the
absorbing set. It is of interest to compute the probability of being absorbed by
a given absorbing set when starting from an initial state i ∈ T . For this we use
the idea of the fundamental matrix, which is related to the number of visits to a
particular state j when starting from a state i.
Theorem 11.2. For a homogenous Markov Chain with transition probability matrix
P , the probability of absorption by the absorbing set starting from transient state i
is
PF = HU (11.3)
170
time to absorption, or time to failure, is taken from the definition of the expectation
of the discrete random variable T ; E(T ), which is given by:
∞
X
Te = E[TiF ] = nP (TiF = n) (11.4)
0
Te = Hc (11.5)
Our task is to model the marker history and lifetime for each device by a
stochastic process with a countable number of states represented in a multistate
model. In general, the future state transitions of a multi-state model may depend
in a complicated way on past events. However, for the special case of a Markov
chain the past and future are independent given its present state. Therefore the
future transitions of a Markov chain depend only on its present state as described
by the transition probabilities Pij (s, t) = P (Y (t) = j|X(s) = i); s < t. We will show
that an estimator for the transition probability matrix can accommodate historical
information, and overcome the apparent limitation inherent to Markov chain models.
Corresponding to the hazard rate for survival, we may for a Markov chain
define the transition intensities
where Y (t− ) denotes the value of the marker Y ”just before” time t. Note that
αij (t)dt is the probability that an device that has experienced i intermittent faults
171
(so is in state i) ”just before” time t, will make a transition to state j in the small
time interval [t, t + dt).
Only for simple Markov chains, is it possible to give explicit expressions for the
probability transition matrix in terms of the transition intensities. For example in
the case of a two-state model, the components of the probability transition matrix
are simply the KM survival and failure estimates. More generally we can express
the (k + 1) × (k + 1) transition probability matrix P (s, t) = {Pij (s, t)} in terms of
the matrix of transition intensities. To see how this is done, we partition the time
interval (s, t] into a number of time intervals s = t0 < t1 < t2 < . . . < tk = t and
use the Markov property to write the transition probability matrix as the product:
If the number of time points increases, while the distance between them goes
to zero uniformly, the matrix product approaches a limit termed a (matrix-valued)
product-integral. The product-integral is written in terms of the (k + 1) × (k +
1) matrix α(u) of transition intensities, that is the matrix where the off-diagonal
elements equal the transition intensities αhj (u), h = j, and the diagonal elements
X
αii (u) = − αij (u)
j=i
are chosen so that all row sums are zero. Since the transition intensities describe the
instantaneous probabilities of transitions between states, P (u, u + du)I + α(u)du,
where I is the (k + 1) × (k + 1) identity matrix. This explains why we may write
Q
the limit of equation (11.7) in product-integral form as: P (s, t) = (s,t] I + α(u)du.
Alternatively, if we let A(t) denote the cumulative transition intensity matrix with
172
Rt
elements Aij (t) = 0
αij (u)du we have:
Y
P (s, t) = {I + dA(u)} (11.8)
(s,t]
P
In the discrete case, the cumulative transition intensity takes the form u≤t
where dNij (t) = Nij (t + dt) − Nij (t), is the increment in the number of transitions
from state i to state j observed over a small time interval [t, t + dt), and R(t) =
{#devices : Ti > t} is the risk set; the number of surviving devices at time t.
P
Furthermore we introduce Âii (t) = − j6=i Âij . The relation in equation (11.8)
suggests that we estimate the matrix of transition probabilities by the k × k matrix
P (s, t), called the Aalen-Johansen (AJ) estimator, Â = {Âij }
Y
P̂ (s, t) = (I + dÂ(u)) (11.10)
(s,t]
The NA estimators are step-functions with a finite number of jumps on (s, t].
Therefore, the AJ estimator given in (11.10) is a finite product of matrices. If one or
more transitions are observed at any time u, then the contribution to (11.10) from
this time point is a matrix I + ∆Â(u), where ∆Â(u) is the k × k matrix with entry
(i, j) equal to ∆Nij (u)/Ri (u) for h 6= j and entry (i, i) equal to −∆Nio (u)/Ri (u)
173
P
with Nio = j6=i Nij . The transition probability matrix defines the instantaneous
probability that an intermittent fault or failure may occur and is estimated using
observations on the marker variable, in this case the event of an error message.
In this section we present a ”worked out” example for evaluating the probabil-
ity transition matrix. We assume that the life-history of the system is described by
a Markov process with a finite number of states Ω = 1, . . . , K. The transition prob-
ability is denoted, as before, by pij (s, t) and describes the probability that a system
in state i at time s transitions into state j by time t. The K × K matrix P (s, t)
summarizes the transition probabilities of the Markov process. We next define the
transition probability matrix in terms of the AJ estimator.
The AJ estimator for the transition matrix P (s, t) is given in equation (11.10).
The interpretation of equation (11.10) is given by equation (11.11), which is the
product over transition matrices between times s and t.
Y
P̂ (s, t) = (I + α̂k ) (11.11)
k:s<tk ≤t
Here,α̂k is a K ×K matrix with entry (i, j) equal to α̂ijk = dijk /Rik , entry (i, i)
equal to α̂iik = −dik /Rik , and all other entries are zero. I is the indentity matrix.
174
The product is taken over increasing tk ’s. In the case of a two-state model, where
state 1 represents the healthy state and state 2 the failure state, the transition and
intensity matrices are two dimensional. In this case the AJ estimator can be seen as
a matrix version of the KM estimator: Ŝ(t) = P̂ [1, 1]. The square brackets indicate
the AJ matrix as opposed to the parenthesis which denote the time components.
Y
P̂ (0, t) = (I + α̂k ) =
k:tk ≤t
dik dijk
Y 1 0 −
= + Rik Rik =
k:tk ≤t 0 1 0 0
175
Tab. 11.1: Lifetime data and survival probability estimates
ID tk dijk From To Rik Ŝ(tk ) F̂ (tk )
13 0.5 0 1 1 16 1 0
1 0.75 1 1 2 15 0.933 0.066
14 0.8 0 1 1 14 0.933 0.066
2 0.91 1 1 2 13 0.861 0.138
3 1.32 1 1 2 12 0.789 0.21
4 1.7 1 1 2 11 0.717 0.282
15 1.7 0 1 1 11 0.652 0.347
16 2.08 0 1 1 9 0.652 0.347
5 2.15 1 1 2 8 0.571 0.428
6 2.76 1 1 2 7 0.489 0.51
7 2.88 1 1 2 6 0.407 0.592
8 2.98 1 1 2 5 0.326 0.673
9 4.51 1 1 2 4 0.244 0.755
10 6.23 1 1 2 3 0.163 0.836
11 8.57 1 1 2 2 0.081 0.918
12 10.23 1 1 2 1 0 1
dik dijk
Y 1−
= Rik Rik =
k:tk ≤t 0 1
Q
dik Q dijk
1− 1− 1−
= k:tk ≤t
Rik k:tk ≤t Rik =
0 0
Ŝ(t) 1 − Ŝ(t)
=
0 1
Using the above data one can show that the following results hold
14 1
15 15 0.933 0.066
P̂ (0, t1 ) = (I + α̂1 ) = =
0 1 0 1
14 1 12 1
15 15 13 13
P̂ (0, t2 ) = P̂ (0, t1 )(I + α̂2 ) = =
0 1 0 1
176
56 9
65 65 0.861 0.138
= =
0 1 0 1
56 9 11 1
65 65 12 12
P̂ (0, t3 ) = P̂ (0, t2 )(I + α̂3 ) = =
0 1 0 1
154 41
195 195 0.789 0.21
= =
0 1 0 1
The AJ survival probability estimates are plotted in Fig. 11.4. The survival
probability is based on the (1,1) element of the AJ estimator at each time step tk ,
as seen from the results above. The resulting survival probability plot in figure
(Fig. 11.4) is a step function which drops at each observed failure time. From figure
11.4 we can infer the survival probability of a new device surviving at time t∗ . For
example, a new device, surviving at t∗ = 4, has a 65% chance of failing in the
next time instance. To summarize, the AJ estimator, non-parametrically estimates
the survival probability as a function of time by taking a product of the transition
intensities of a Markov chain at each of the observed event times. Because here we
only have two types of events, there is only one transition: from healthy to failure.
To see the usefulness of this model, we next present a case study that generalizes
these computations to higher dimensions to accommodate recurrent intermittent
faults.
The expected time to failure is computed using Markov chain theory on ab-
sorbing chains as discussed in section 12.4. Specifically, we are interested in the
fundamental matrix of the transition probability matrix H. As stated in 11.3, the
ij entry of H gives the expected number of steps the chain takes to reach state j
if it starts in the state i. Consider for example, the transition probability matrix
given by the AJ estimator at time t3 , in this case the transient state matix Q is
177
1
P(0,t1) AJ Estimator
P(0,t2)
0.8
Survival Probability
0.6
0.4
0.2
0
0 2 4 6 8 10 12
Time
The first entry in the expected time to failure vector is the survival estimate
at the earliest lifetime measurement t1 = 0.5, which is a censoring time. Because at
this time, no failures have previously been observed, this censored lifetime does not
contribute to the survival probability estimate. The AJ estimate at time t3 , is the
survival probability estimate at the third failure time, or including censored times,
its the fifth lifetime measurement, as seen in the vector above.
In a higher dimensional multistate model, we need to take the inverse of Q
178
14
13
12
11
10
matrices, which becomes a computational problem when there are many states to
account for. Difficulties may arise in taking the inverse of anon-singular matrix. Here
we assume that computational techniques for matrix inverse evaluation are available
(LU decomposition, Singular Value decomposition, Gauss-Jordan elimination, etc.).
In the following case study the same calculations are used but in matrix form.
In tests, failures are repeatedly induced by stressing the computer with sim-
ulated usage. Failures are defined as unwanted/automatic restarts or shutdowns.
Computer usage is simulated by running intensive programs that consume compu-
tational resources, and cause the computer to freeze or trigger a shutdown. It is
believed that errors are correlated to failure and can therefore be used as precur-
sors/markers in our multistate model. Table 11.2 shows a sample of failure times
observed for a computer under test. The first column shows the start time after a
restart/reboot of the computer. The second column is the time the computer was
observed to fail: time to failure (TTF). The third column records the calendar time
179
Tab. 11.2: Sample of failure times for collected from a computer
Start Stop TTF TBF
2/23/2010 7:53 2/24/10 11:05 4:33 27.19
2/24/2010 11:05 2/25/10 7:26 12:57 20.35
2/25/2010 7:26 2/25/10 8:09 6:14 0.73
2/25/2010 8:09 2/26/10 13:28 13:55 29.32
2/26/2010 13:28 2/26/10 23:16 9:07 9.79
in hours from the start of the experiment (t=0). The last column shows the time in
hours between each failure (TBF). This is the lifetime of an device as discussed in
section 3.
In addition to failure time information, we also collect error event times. This
data is collected in a similar way. Lifetimes and error event times are collected into a
”life table”, which contains all the data: lifetimes and error event times in one table
structure. Table 11.3 gives an example of a life table. Each row in the life table
represents the occurrence of an event: failure or error. In Table 11.3, the failure
ID represents an device in the computer. There may be multiple error events prior
to each failure; therefore, the same failure ID can appear multiple times. When a
failure is observed, the failure indicator If = 1, otherwise If = 0 to indicate an error
event.
We simulated lifetime and error data. The main objective of this case study is
to test and validate the proposed multistate model. Our simulation design hypothe-
sis is: 1) that with each error, the probabilities of another error or failure occurring
increase and 2) with each failure the probability that an error occurs increases. In
180
other words, devices that experience more errors typically fail earlier, and a com-
puter which has seen many failures is more likely to generate more errors than a
computer which has not seen as many failures.
In our simulation, random lifetimes and error times are generated based on
the FHT model discussed in section 12.4. Lifetimes Tq , therefore, are modeled as
the first hitting time of a degradation process X(t) to a fixed failure threshold a,
namely, T = inf t : X(t) = a. We modeled the degradation process by a Wiener
process with drift given by:
181
1
0.8
0.6 Degradation
X(t),Y(t)
Process X(t)
0.4
Event/Marker
0.2 Process Y(t)
0
0 5 10 15 20 25 30 35 40
Time
Fig. 11.6: Simulation of the degradation process conditioned on an error event process
182
350
No Errors
300
250
Frequency
200
150
100
50
0
20 40 60 80 100 120 140 160 180 200
Expected Time to Failure
250
With 1 Error
225
200
175
150
Frequency
125
100
75
50
25
0
40 45 50 55 60 65 70 75 80 85 90 95 100
Expected Time to Failure
450
Two Errors
400
350
300
250
200
150
100
50
0
20 30 40 50 60 70 80 90 100 110 120 130 140 150 160
Expected Time to Failure
Fig. 11.7: Expected time to failure distribution for devices with 0,1 and 2 error events
for devices that experience 2 errors. Using random lifetime and error times gener-
ated from our statistical model this result validates the simulation design discussed
earlier.
183
p01 p01 0 p03
0 p11 p12 p13
0 0 p22 p23
0 0 0 1
The multistate model, in this case, consists of two marker states and a failed
state. In this model only two errors occurrences are considered prior to failure
for each device, but this need not be the case, any discrete number of errors can
be accommodated. State 0 represents the state in which the computer has not
experienced any errors, state 1 the state in which the computer experiences its first
error and state 2 its second and last error. States 0 through 2 represent the transient
states in a Markov chain. State 3, the failure state, represents the absorbing state.
Fig. 11.8 shows the structure of the transition probability matrix used in this case
study for a three state multistate model.
In a three state multistate model, each device under observation has a maxi-
mum of three possible event times, two for the errors and one for the devices failure.
For example, ID=4 experiences its first error at 29.89 hours, its second at 32.82
hours and fails 41.37 hours. Time is measured from t=0 for all devices. For ID=1,
it does not experience any errors and it fails at 59.63 hours.
The time considered in this experiment starts at t= 0 and ends at the time
of the last observed failure across the 40 lifetimes, which for this data set is t=
83.36 hours. The AJ estimator is applied over a time increment of r= 0.00018 hours
which is equal to the smallest time difference between any two error events in the
data set. This guarantees that at any time, at most, only 1 error event can occur.
The transition probability matrix therefore is estimated at every time step, in total
4631 times, to generate a three dimensional matrix k × k × m, with k = 4 and
m = 4631. Each slice of this matrix along m, gives the non-parametric estimate of
the transition probability matrix at a particular time.
184
1
0.9 From State 0
0.8 From State 1
Survival Probability
0.7 From State 2
0.6
0.5
0.4
0.3
0.2
0.1
0
0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 84
Time (Hours)
Fig. 11.9: Aalen-Johansen survival probability estimate starting from states 0, 1 and 2
Figure 11.9 shows the AJ estimator results for the survival probabilities start-
ing from states 0, 1 and 2. These plots are generated based on a Matlab code
developed at CALCE to implement the AJ estimator (equation (11.10)). The eval-
uation of the AJ estimator is performed at discrete times, similar to example 5.1.
Survival probability from state 0 is based on element (1,4) in the transition proba-
bility matrix (see figure 11.6), from state 1, based on element (2,4), and from state
2 based on element (3,4).
From the results in figure 11.9 we observe that the AJ estimator detects the
design conditions in the simulated data. We can see that devices that experience
an increasing number of errors before failure have increasingly lower chances of
surviving, and this result validates the simulation design discussed earlier.
Therefore, for a new computer, if we know how long it’s been running since
last shutdown, and we know the number or errors in any (errors of some predeter-
mined type) it experienced, we can get point estimates of its remaining life. For
example, for a test computer that is surviving at time t, its expected time to failure
is plotted in Fig. 11.10, for times t ranging from 0 to 84 hours. Using equation
(11.5), figure 11.10 plots the expected failure-time as the first hitting time of the
185
80 From State 0
From State 1
70 From State 2
50
40
30
20
10
0
15 20 25 30 35 40 45 50 55 60 65 70 75 80 85
Time (Hours)
Fig. 11.10: Expected time to failure starting from state 0,1 and 2 respectively
186
a list of all possible points in the bivariate state space. Each state represents the
joint degradation/marker values indexed by time. The marker value is observable
and represents the cumulative number of errors experienced by the device, while
the degradation variable is unobservable and therefore unknown. In this approach,
inference is based on the progress of the marker variable which is assumed correlated
to the degradation variable.
Inference of the transition probability matrix is approached non-parametrically
using the AJ estimator over small time increments spanning the life time of the
device. A reducible absorbing time homogeneous Markov chain is used to compute
the expected first hitting times to estimate the expected time to failure. A case
study simulates a sample of 40 lifetimes and 45 error event times. The simulation
generates failure times and event times according to a design hypothesis discussed
in the chapter, and used to emulate experimental data. In the simulation, the
marker variable is modeled as a Poisson random variable, and the degradation as
a Gaussian random variable. The degradation process is conditioned on the event
process; its drift parameter increases with increasing number of errors and failures in
the computer. From the results we see, as anticipated, the expected time to failure
decreases given an increasing number of errors experienced by an device.
This chapter contributes to the literature on degradation models on the fol-
lowing levels:
187
process are summarized using a multistate representation which is simple and
intuitive to use for prediction.
3. Using counting process theory, time dependent covariates (in this case the
occurrence of an error) influence the estimation of transition probabilities.
This is perhaps the most useful aspect of the model because it can naturally
accommodate a large number of covariates (in this thesis we use one) without
burdening computations and complicating predictive inference equations.
5. Connects the model with PHM, a methodology that requires real-time health
assessment and predictions. Extensions to this model can include a parametric
model to describe the relationship between the marker and a degradation vari-
able. In this case, to achieve a valuable model, the degradation variable must
represent the mode of a known failure mechanism, and the marker variable
should be correlated to the degradation variable.
188
Tab. 11.4: Simulated Data for Case Study
In order to predict the remaining life of a device we need to know when other
similar devices failed, as well as how they responded to stress over time. To get
this information we must conduct tests, and expose a sample of devices to stress
and measure their responses. Only then can we take information from a new fielded
device and make any inference on its reliability. In this work the au The impetus for
degradation models in PHM stems from the need to explain heterogeneous reliability
qualities across a sample of devices used in dynamic stress environments. This need
is further exacerbated by the requirement in PHM to predict failure-times when
failure-time samples are small, and when degradation data are not predictive of fail-
ure. Small failure-time samples are common in highly reliable products or products
with short product-cycles. Small failure-time samples also result from reduced ac-
celerated test conditions, aimed to preserve the failure generating mechanisms for
devices put through failure-tests.
Degradation data are collected for each device and used as auxiliary reliability
information to improve reliability models. Heterogeneous degradation data, col-
lected from degradation variables, can provide valuable insight into device-specific
reliability. First hitting time models use degradation variables to define the failure-
time, implicitly enforcing a causality between the underlying failure mechanism,
degradation and the failure-time. Degradation variables play another important
role in enhancing reliability and PHM models, because, they represent responses to
stress or usage, which varies across devices. In this way degradation variables allow
us to model the effect of dynamic environments on our failure-time predictions.
Central to the thesis are so called non-predictive degradation variables. These
are known, typically observable, degradation variables that attain the failure -
threshold level suddenly without any preceding trend, a trend useful for making
predictions on lifetimes. In such a case, as we discuss in this thesis, we use latent
degradation models with terminal observations on an otherwise latent true degra-
dation variable and longitudinal measurements on observable marker variables. We
provide justification for using bivariate latent degradation models for data collected
in failure-tests, and through our model we address key data-limitations encountered
in such settings.
Our baseline analytical framework is Whitmore’s bivariate Wiener model for
terminal degradation and marker data-observations. In chapter 1 we introduce the
problem and motivate a direction of research. In chapter 2 we present the data-
structures used in the thesis and examples taken from failure-tests of electronics
conducted at CALCE. In chapter 3 we present main theorems and lemmas on es-
timation theory. In chapters 4, 5 and 8 we present extensions to Whitmore’s first
hitting time model. In chapter 6 we present case studies that compare the perfor-
mance of our degradation model on various data-structures, the results of which
form the basis of our contributions. In chapter 7 we consider the effect of covariates
on estimation under the same degradation model.
The first extension and first contribution of this work is to address and provide
a simple solution to the ”small failure-time sample” problem. We developed para-
metric and predictive inference equations based on Whitmore’s first hitting time
model, for a data-structure that augments terminal degradation measurements to
a terminal data-structure. With terminal degradation observations on both failed
and surviving devices we were able to use, in contrast to Whitmore, terminal degra-
dation information on surviving devices. In other words, in our approach, surviving
devices become much more valuable for estimation and inference.
191
We compared the mean failure-time parameter of an IG lifetime distribution,
and we showed that our approach under the TMD data-structure consistently re-
duces the asymptotic variance of MLEs. More importantly, our experimental results
indicate that our approach performs increasingly better under smaller failure-time
sample sizes. As expected, from a statistical point of view, the efficacy of the TMD
over the TM data-structure is explained because ML estimation improves with en-
hanced data, when sample sizes are large.
The improvement in inference under the TMD data-structure also has broader
commercial implications. The results show that we can ”afford” to reduce designed
accelerating conditions, and test durations in planned failure-tests. From an engi-
neering perspective, reducing accelerating conditions is a welcome option, because,
as we discuss in the introduction, highly accelerated test conditions may alter tar-
geted failure-mechanisms. Typically, failure-tests are conducted over a time-period
long enough to see enough devices fail. Shorter test-durations are not only less
costly, they are sometimes required due to short product life-cycles. Under latent
degradation conditions, our results suggest investing in failure-analysis equipment
capable of measuring terminal degradation.
The second contribution of the thesis lies in our treatment of longitudinal
marker observations in a latent degradation model for both parametric and predic-
tive inference. We develop parametric inference for a general longitudinal marker
data-structure and show in our results improved estimation starting from two in-
termediate marker observations. The efficacy of the longitudinal data-structure
depends on the strength of correlation between the marker and the degradation
variable. In our simulations and analysis, we use a strong correlation coefficient.
Further work is needed to test estimation improvement under weaker correlations.
The third contribution of the thesis incorporates a variable failure-threshold
to Whitmore’s bivariate Wiener model. We develop parametric inference equations
192
with the failure-threshold variable modeled as a Gaussian random variable and in-
dependent of the degradation process. More realistically, as part of future work,
the failure-threshold variable may also be dependent upon the degradation process,
and can itself be modeled, more generally as a stochastic process. We are currently
working on code to evaluate AREs under the TMD data-structure for fixed vs vari-
able threshold models. These results are not availabe in the thesis, but are however
anticipated to be part of subsequent publications.
In chapter 7, in addition to presenting covariates as part of FHT models, we
present machine learning approaches as part of degradation models. Specifically,
we consider Gaussian process regression and support vector machine classification
as methods for including covariates. Our introduction, and justification of SVMs in
this context constitutes our fourth and last contribution in part-I of the thesis.
In part-II of the thesis we investigate degradation models based on nonpara-
metric models. In chapter 9 we analyze the reliability of a product from a health
monitoring perspective in the context of PHM. In the absence of failure training-
data, anomaly detection is approached through a one-class learning algorithm based
on SVM classification. This is also used in a Bayesian framework to estimate the
posterior class probabilities of test data with unknown class. In this work we make
a contribution to the field of reliability by interpreting health as the outcome of a
classifier. We introduce a methodology that connects machine-learning analysis to
FHT models by helping determine suitable marker variables that can be used to
track latent degradation. We solve a novelty detection problem with a one-class
classification algorithm, and a Bayesian framework for uncertainty analysis. The
results of our Bayesian classifier, we argue and show through case studies, are suited
for further trending and event-time predictions.
In chapter 11 we consider degradation variables that take-on discrete values,
and use computer reliability as a working example. Specifically we are interested
193
in modeling the occurrence of intermittent errors that lead to failure. Here we use
a multistate Markov model to develop a methodology to model failure-time data
together with event-time data (errors). Failure-times are modeled as the FHT to
a fixed failure state in a multistate model. The marker value is observable and
represents the cumulative number of error experienced by the device, while the
degradation variable is unobservable and therefore unknown. In this approach,
inference is based on the process of the marker variable which is assumed correlated
to the degradation variable.
Inference on the transition probability matrix is approached using the non-
parametric Aalen-Johansen estimator over small time increments spanning the life-
time of the device. A reducible absorbing time-homogeneous Markov chain is used to
compute the expected first hitting times. A simulation case-study is used to generate
failure-times an event times. From the results we observe, that as anticipated, the
expected failure-time decreases given an increasing number of errors experienced by
a device. This part of the thesis contributes to the literature on degradation models
on the following levels: (i) it casts a degradation model as a multistate model, (ii) it
accommodates covariates using counting process theory, (iii) it avoids making many
parametric assumptions and (iv) it connects the model with PHM.
194
BIBLIOGRAPHY
[1] Albert, P.S., Longitudinal data analysis (repeated measures) in clinical trials,
Statistics in Medicine, Vol. 18, No. 13, 1707–1732, 1999
[4] M. Pecht, A. Dasgupta, D. Barker, C.T. Leonard, The reliability physics ap-
proach to failure prediction modelling, Quality and Reliability Engineering In-
ternational, Vol. 6, No. 4, 267–273, 1990
[5] J.M. Hu, M. Pecht, A. Dasgupta,A probabilistic approach for predicting ther-
mal fatigue life of wire bonding in microelectronics, Journal of Electronic Pack-
aging, Vol. 113, 275–285, 1991
[10] Dasgupta A., M. Pecht, Material failure mechanisms and damage models, IEEE
Transactions on Reliability, Vol. 40, No. 5, 531–536, 1991
[14] J. Lu, Degradation processes and related reliability models, Thesis, 1995
196
[17] M. Pecht, A prognostics and health management roadmap for information and
electronics-rich systems, Microelectronics Reliability, Vol. 50, No. 3, 317–323,
2010
[19] W. Kahle, Simultaneous confidence regions for the parameters of damage pro-
cesses, Statistical Papers, Vol. 35, No. 1, 27–41, 1994
[21] J. Lu, Degradation processes and related reliability models, Unpublished thesis,
McGill University, 1995
[23] K.A. Doksum, S.L.T. Normand, Gaussian models for degradation processes-
Part I: Methods for the analysis of biomarker data, Lifetime Data Analysis,
Vol. 1, No. 2, 131–144, 1995
[25] L.I. Pettit, K.D.S. Young, Bayesian analysis for inverse Gaussian lifetime data
with measures of degradation, Journal of statistical computation and simula-
tion, Vol. 63, No. 3, 217–234, 1999
197
[26] M. Lee, V. DeGruttola, D. Schoenfeld, A model for markers and latent health
status, Journal of the Royal Statistical Society: Series B (Statistical Method-
ology), Vol. 62, No. 4, 747–762, 2000
[28] M.L.T. Lee, G.A. Whitmore, Threshold regression for survival analysis: model-
ing event times by a stochastic process reaching a boundary, Statistical Science,
Vol. 21, No. 4, 501–513, 2006
[29] J. Tang, T.S. Su Estimating failure time distribution and its parameters based
on intermediate data from a Wiener degradation model, Naval Research Logis-
tics, Vol. 55, No. 3, 265–276, 2008
[31] [Link], W. Wang, C.H. Hu, D.H. Zhou, Remaining useful life estimation-A re-
view on the statistical data driven approaches, European Journal of Operational
Research, Vol. 213, No. 1,1–14, 2010
[33] M.L.T. Lee, G.A. Whitmore, B.A. Rosner, Threshold regression for survival
198
data with time-varying covariates, Statistics in medicine, Vol. 29, No. 7-8, 896–
905, 2010
[36] N.D. Singpurwalla, Inference from accelerated life tests when observations are
obtained from censored samples, Technometrics, Vol. 13, No. 1, 161–170,1971
[38] N.D. Singpurwalla Inference from Accelerated Life Tests Using Arrhenius Type
Re-Parameterizations, Technometrics, Vol. 15, No. 2, 289–299, 1973
[39] N.R. Mann, R.E. Schafer, N.D. Singpurwalla, Methods for statistical analysis
of reliability and life data, Wiley, 1974
[40] G.K. Bhattacharyya, A. Fries, Inverse Gaussian regression and accelerated life
tests, Lecture Notes-Monograph Series, Vol. 2, 101–117, 1982
[41] W. Nelson, Accelerated testing: statistical models, test plans, and data analy-
ses, Wiley, 1990
[42] M.B. Carey, R.H Koenig, Reliability assessment based on accelerated degrada-
tion: a case study, IEEE Transactions on Reliability, Vol. 40, No. 5, 499–506,
1991
199
[43] K.A. Doksum, A. Hóyland, Models for variable-stress accelerated life testing
experiments based on Wiener processes and the inverse Gaussian distribution,
Technometrics, Vol. 34, No. 1, 74–82, 1992
[44] W.Q. Meeker, L.A. Escobar, A review of recent research and current issues in
accelerated testing, International Statistical Review/Revue Internationale de
Statistique, Vol. 61, No. 1, 147–168, 1993
[47] W.Q. Meeker, L.A. Escobar, C.J. Lu, Accelerated degradation tests: modeling
and analysis, Technometrics, Vol. 40, No. 2, 89–99, 1998
[48] W.J. Owen, W.J. Padgett, Accelerated test models for system strength based
on Birnbaum-Saunders distributions,Lifetime Data Analysis, Vol. 5, No. 2, 133–
147, 1999
[49] A. Onar, W.J. Padgett, Accelerated test models with the inverse Gaussian
distribution, Journal of statistical planning and inference, Vol. 89, No. 1-2,
119–133, 2000
200
[51] V. Bagdonavicius, O. Cheminade, M. Nikulin, Statistical planning and inference
in accelerated life testing using the CHSS model, Journal of statistical planning
and inference, Vol. 126, No. 2, 535–551, 2004
[52] W.J. Padgett, M.A. Tomlinson, Inference from accelerated degradation and
failure data based on Gaussian process models, Lifetime Data Analysis, Vol.
10, No. 2, 191–206, 2004
[53] C. Park, W.J. Padgett, Accelerated degradation models for failure based on
geometric Brownian motion and gamma processes, Lifetime Data Analysis,
Vol. 11, No. 4, 511-527, 2005
[55] S.J. Bae, W. Kuo, P.H. Kvam, Degradation models and implied lifetime dis-
tributions,Reliability Engineering & System Safety, Vol. 92, No. 5, 601–608,
2007
[56] W.Q. Meeker, L.A. Escobar, Y. Hong, Using accelerated life tests results to
predict product field reliability, Technometrics, Vol. 51, No. 2, 146–161, 2009
201
[59] S. Han, M. Osterman, M. Pecht, Electrical Shorting Propensity of Tin
Whiskers, IEEE Transactions on Electronics Packaging Manufacturing, Vol.
33, No. 3, 205–211
[62] W.H. Press, B.P. Flannery, [Link], W.T. Vetterling, Numerical Recipes,
Cambridge University Press, 1992
[63] , L.B. Koralov, Y.G. Sinai, Theory of probability and random processes,
Springer Verlag, 2007
[65] M.L. Eaton, Multivariate statistics: a vector space approach, Wiley New York,
1983
[68] , R.S. Chhikara, J.L. Folks, The inverse Gaussian distribution as a lifetime
model,Technometrics, Vol. 19, No 4. 461–468, 1977
202
Transactions on Components and Packaging Technologies, Vol. 29, No. 1, 1521–
3331, 2006
[72] X.K. Song, Correlated data analysis: modeling, analytics, and applications,
Springer, 2007
[73] N.M. Laird, J.H. Ware, Random-effects models for longitudinal data, Biomet-
rics Vol.38, 963–974, 1982
[74] J.P. Klein, M.L. Moeschberger, Survival analysis: techniques for censored and
truncated data, Springer Verlag, 2003
203
ogy for electronic systems”, IEEE Transactions on Components and Packaging
technologies, vol. 26, no. 3, pp. 625–634, 2003
[80] N. Vichare, P. Rodgers, and V. Eveloy, and M. Pecht, ”In Situ Temperature
Measurement of a Notebook ComputerA Case Study in Health and Usage Mon-
itoring of Electronics”, IEEE Transactions on Device and Materials Reliabil-
ity,vol. 4, no. 4, pp. 658–663, 2004
[81] T. Stibor, P. Mohr, and J. Timmis, and C. Eckert, ”Is negative selection ap-
propriate for anomaly detection?”, in Proceedings of the 2005 conference on
Genetic and evolutionary computation, ACM, pp. 321–328, 2005
[83] D. Tax and P. Juszczak, ”Kernel whitening for one-class classification”, Inter-
national Journal of Pattern Recognition and Artificial Intelligence, vol. 17, no.
3, pp. 333–347, 2003
[84] H. Wang, and Z. Song, and P. Li, ”Fault detection behavior and performance
analysis of principal component analysis based process monitoring methods”,
Ind. Eng. Chem. Res, vol. 41, no. 10, pp. 2455–2464, 2002
204
[85] J. Jackson and G. Mudholkar, ”Control procedures for residuals associated with
principal component analysis”, Technometrics, vol. 21, no. 3, pp. 341–349, 1979
[86] V. Klema and A. Laub, ”The singular value decomposition: Its computation
and some applications”, IEEE Transactions on Automatic Control, vol. 25, no.
2, pp. 164–176, 1980
[87] L. Ruixin, W. Dongfeng, and H. Pu et al., ”On the applications of SVD in fault
diagnosis”, in IEEE International Conference on Systems, Man and Cybernet-
ics, vol 4, no. 5, pp. 3763–3768, 2003
[88] V. Vapnik, ”The nature of statistical learning theory”, Springer Verlag, 2000
[93] C. Hsu, C. Chang, and C. Lin et al., ”A practical guide to support vector
classification”, Citeseer, 2000
205
of SVMs with an application to unbalanced classification”,Advances in Neural
Information Processing Systems, vol. 18, pp. 467–474, 2006
[95] J. Kwok, ”The evidence framework applied to support vector machines”, IEEE
Transactions on Neural Networks, vol. 11, no. 5, pp. 1162–1173, 2000
[97] J. Kwok, ”Moderating the outputs of support vector machine classifiers”, IEEE
Transactions on Neural Networks, vol. 10, no. 5, pp. 1018–1031, 1999
[99] W. Fan, M. Miller, S. Stolfo, W. Lee, W. and P. Chan, ”Using artificial anoma-
lies to detect unknown and known network intrusions”,Knowledge and Infor-
mation Systems, vol. 6, no. 5, pp. 507–527, 2004
[101] H. Yu, J. Han, and K. Chang, ”PEBL: Web page classification without neg-
ative examples”,IEEE Transactions on Knowledge and Data Engineering, vol.
16, no. 1, pp. 70–81, 2004
[102] J. Theiler and D. Cai, ”Resampling approach for anomaly detection in multi-
spectral images”,in Proceedings of SPIE, vol. 5093, pp. 230–240, 2003
206
[103] S. Marsland, ”Novelty detection in learning systems”, Neural computing sur-
veys, vol. 3, pp. 157–195, 2003
[109] J. Nuñez Garcia, Z. Kutalik, K. Cho, and O. Wolkenhauer, ”Level sets and
minimum volume sets of probability density functions”, International journal
of approximate reasoning, vol. 34, no. 1, pp. 25–47, 2003
207
kernel for SVM classification in multimedia applications”,Advances in Neural
Information Processing Systems, vol. 16, 2004
[112] J. Gao and P. Tan, ”Converting output scores from outlier detection algo-
rithms into probability estimates”, in ICDM06: Proceedings of the Sixth Inter-
national Conference on Data Mining, pp. 212–221
[113] J. Platt, ”Probabilistic outputs for support vector machines and comparisons
to regularized likelihood methods”, Advances in large margin classifiers, pp.
61–74, 1999
[114] W. Chu, S. Keerthi, and C. Ong, ”A new Bayesian design method for support
vector classification”, in Special Section on Support Vector Machines of the 9th
International Conf. on Neural Information Processing, Citeseer, pp. 888–892,
2002
[115] Sotiris, V.A. and Tse, P.W. and Pecht, M.G., ”Anomaly Detection Through
a Bayesian Support Vector Machine”,Reliability, IEEE Transactions on, Vol.
59,no. 2,pp. 277–286,2010
[116] Doksum, K.A. and Hóyland, A., ”Models for variable-stress accelerated life
testing experiments based on Wiener processes and the inverse Gaussian dis-
tribution”, Technometrics,no. 1, 74–82,1992
[117] Doksum, K.A. and Normand, S.L.T., ”Gaussian models for degradation
processes-Part I: Methods for the analysis of biomarker data”,Lifetime Data
Analysis, Vol. 1,no. 2, pp. 131–144,1995
208
[119] Satten, G.A. and Longini Jr, I.M., ”Markov Chains With Measurement Error:
Estimating theTrue’Course of a Marker of the Progression of Human Immun-
odeficiency Virus Disease”, Applied Statistics, vol. 45,no. 3,pp. 275–309, 1996
[120] Henderson, R. and Diggle, P. and Dobson, A., ”Joint modelling of longitudinal
measurements and event time data”, Biostatistics, vol. 1, no. 4,pp. 465–480,
2000
[121] Meeker, W.Q. and Escobar, L.A., ”Statistical methods for reliability data”,
Wiley New York, 1998
[124] Bagdonaviius, V. and Bikelis, A. and Kazakeviius, V. and Nikulin, M., ”Non-
parametric estimation from simultaneous degradation and failure time data”,
Comptes Rendus Mathematique, vol. 335, no. 2, pp. 183–188, 2002
[125] Putter, H. and van der Hage, J. and de Bock, G.H. and Elgalta, R. and van de
Velde, C.J.H., ”Estimation and Prediction in a Multi-State Model for Breast
Cancer”,Biometrical journal, vol. 48,no. 3, pp. 366–380, 2006
[126] Peña, E.A., ”Dynamic modelling and statistical analysis of event times”, Sta-
tistical science: a review journal of the Institute of Mathematical Statistics,
vol. 21,no. 4, pp. 487–500, 2006
209
[127] Meira-Machado, L. and de Uña-Álvarez, J. and Cadarso-Suárez, C. and Ander-
sen, P.K., ”Multi-state models for the analysis of time-to-event data”,Statistical
methods in medical research, vol. 18, no. 2, pp. 1–28,2009
[129] Andersen, P.K. and Pohar Perme, M., ”Inference for outcome probabilities in
multi-state models”, Lifetime data analysis, vol. 14, no. 4, pp. 405–431, 2008
[130] Cook, R.J. and Lawless, J.F. and Lakhal-Chaieb, L. and Lee, K.A., ”Ro-
bust estimation of mean functions and treatment effects for recurrent events
under event-dependent censoring and termination: application to skeletal com-
plications in cancer metastatic to bone”, Journal of the American Statistical
Association, vol. 104, no. 485, pp. 60–75, 2009
[131] Aalen, O.O. and Borgan, Ø. and Gjessing, H.K., ”Survival and event history
analysis: a process point of view”, Springer Verlag, 2008
[134] Aalen, O.O. and Johansen, S., ”An Empirical Matrix for Non-Homogeneous
Markov Chains Based on Censored Observations”, Scandinavian Journal of
Statistics, Vol. 5, No. 3, pp. 141–150, 1978
210
[135] Brémaud, P., ”Markov chains: Gibbs fields, Monte Carlo simulation, and
queues”, Springer, 1999
[136] Singpurwalla, N.D., ”On competing risk and degradation processes”, Lecture
Notes-Monograph Series, vol. 49, pp. 229–240, 1996
211