Notes on Statistical Methods for Reliability
Data
Qianqian Shan
Contents
1 Reliability Concepts and Reliability Data 7
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2 Genral Models for Reliability Data . . . . . . . . . . . . . . . 7
1.2.1 Target Population or Process . . . . . . . . . . . . . . 7
1.2.2 Causes of Failure and Degradation Leading to Failure . 8
1.2.3 Environment Effects on Reliability . . . . . . . . . . . 8
1.2.4 Time Scale . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3 Repairable Systems and Nonrepairable Units . . . . . . . . . . 9
1.4 Strategy for Data Collection, Modeling and Analysis . . . . . 9
1.4.1 Planning a Reliability Study . . . . . . . . . . . . . . . 9
1.4.2 Strategy for Data Analysis and Modeling . . . . . . . . 10
2 Models, Censoring, and Likelihood for Failure-Time Data 11
2.1 Models for Continuous Failure-Time Process . . . . . . . . . . 11
2.2 Models for Discrete Data From a Continuous Process . . . . . 14
2.3 Censoring . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.3.1 Censoring Mechanism . . . . . . . . . . . . . . . . . . 14
2.3.2 Assumptions on Censoring Mechanism . . . . . . . . . 15
2.4 Likelihood . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3 Non-parametric Estimation 16
3.1 Estimation from Singly Censored Interval Data . . . . . . . . 16
3.2 Basic Ideas of Statistical Inference . . . . . . . . . . . . . . . . 16
3.3 Confidence Intervals from Complete or Singly Censored Data . 16
3.3.1 Point-wise Binomial-based Confidence Interval for F (t
i
) 16
3.3.2 Pointwise Normal-Approximation Confidence Interval
for F (t
i
) . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.4 Estimation from Multiply Censored Data . . . . . . . . . . . . 17
3.5 Pointwise Confidence Intervals from Multiply Censored Data . 17
1
CONTENTS
3.5.1 Approximate Variance of
ˆ
F (t
i
) . . . . . . . . . . . . . . 17
3.5.2 Greenwood’s Formula . . . . . . . . . . . . . . . . . . . 18
3.5.3 Pointwise Normal-Approximation Confidence Interval
for
ˆ
F (t
i
) . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.6 Estimation from Multiply Censored Data with Exact Failures 18
3.7 Simulataneous Confidence Bands . . . . . . . . . . . . . . . . 19
3.8 Uncertain Censoring Times . . . . . . . . . . . . . . . . . . . 19
4 Loaction-Scale-Based Parametric Distributions 19
4.1 Quantities of Interest in Reliability Applications . . . . . . . . 19
4.2 Location-scale and Log-location-scale Distributions . . . . . . 20
4.3 Exponential Distribution . . . . . . . . . . . . . . . . . . . . . 21
4.4 Normal Distribution . . . . . . . . . . . . . . . . . . . . . . . 21
4.5 Log-normal Distribution . . . . . . . . . . . . . . . . . . . . . 21
4.6 Smallest Extreme Value Distribution . . . . . . . . . . . . . . 21
4.7 Weibull Distribution . . . . . . . . . . . . . . . . . . . . . . . 21
4.8 Largest Extreme Value Distribution . . . . . . . . . . . . . . . 22
4.9 Logistic Distribution . . . . . . . . . . . . . . . . . . . . . . . 22
4.10 Log-logistic Distribution . . . . . . . . . . . . . . . . . . . . . 22
4.11 Generating Pseudo-random Observations from a Specified Dis-
tribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4.11.1 Pseudo-random Observations from Continuous Distri-
butions . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4.11.2 Pseudo-random Observation from Discrete Distributions 22
4.11.3 Efficient Generation of Censored Pseudorandom Samples 23
5 Other Parametric Distributions 23
6 Probability Plotting 24
6.1 Purpose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
6.2 Linearing a CDF . . . . . . . . . . . . . . . . . . . . . . . . . 24
6.3 Probability Plotting Positions . . . . . . . . . . . . . . . . . . 24
6.3.1 Criteria for Choosing Plotting Positions . . . . . . . . 24
6.3.2 Choice of Plotting Positions . . . . . . . . . . . . . . . 25
6.4 Probability Plots with Specified Shape Parameter . . . . . . . 25
7 Parametric Likelihood Fitting Concepts: Exponential Distri-
bution 26
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
7.2 Parametric Likelihood . . . . . . . . . . . . . . . . . . . . . . 26
7.3 Confidence Intervals for θ . . . . . . . . . . . . . . . . . . . . 26
2
CONTENTS
7.3.1 Likelihood Confidence Intervals for θ . . . . . . . . . . 26
7.3.2 Normal-Approximation Confidence Intervals for θ . . . 27
7.4 Confidence Intervals for Function of θ . . . . . . . . . . . . . . 27
7.5 Likelihood for Exact Failure Times . . . . . . . . . . . . . . . 28
7.5.1 Correct Likelihood for Observations Reported as Exact
Failures . . . . . . . . . . . . . . . . . . . . . . . . . . 28
7.5.2 Using Density Approximation for Observations Reported
as Exact Values . . . . . . . . . . . . . . . . . . . . . . 28
7.6 Data Analysis with No Failures . . . . . . . . . . . . . . . . . 28
8 Maximum Likelihood for Log-Location-Scale Distributions 29
8.1 Likelihood . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
8.2 Likelihood Confidence Regions and Intervals . . . . . . . . . . 29
8.2.1 Joint Confidence Regions for µ and σ . . . . . . . . . . 29
8.2.2 Individual Confidence Intervals for µ and σ . . . . . . . 29
8.2.3 Likelihood Confidence Intervals for Functions of µ and σ 29
8.3 Normal Approximation Confidence Intervals . . . . . . . . . . 30
8.4 Estimation with Given σ . . . . . . . . . . . . . . . . . . . . . 30
9 Bootstrap Confidence Intervals 30
9.1 Bootstrap Sampling . . . . . . . . . . . . . . . . . . . . . . . . 31
9.2 Confidence Intervals . . . . . . . . . . . . . . . . . . . . . . . 31
9.3 Percentile Bootstrap Method . . . . . . . . . . . . . . . . . . . 32
10 Planning Life Test 33
10.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
10.2 Approximate Variance of ML Estimators . . . . . . . . . . . . 34
10.2.1 Basic Large-Sample Approximation . . . . . . . . . . . 34
10.3 Sample Size for Unrestricted Functions . . . . . . . . . . . . . 35
11 Parametric Maximum Likelihood: Other Models 35
11.1 Truncated Data and Truncated Distributions . . . . . . . . . . 35
11.2 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
12 Prediction of Future Random Quantiles 36
12.1 Probability Prediction Intervals(θ Given) . . . . . . . . . . . . 36
12.2 Statistical Prediction Interval(θ Estimated) . . . . . . . . . . 37
12.2.1 Coverage Probability Concepts . . . . . . . . . . . . . 37
12.2.2 Naive Method for Computing a Statistical Prediction
Interval . . . . . . . . . . . . . . . . . . . . . . . . . . 37
12.3 The (Approximate) Pivotal Method for Prediction Intervals . 38
3
CONTENTS
12.3.1 Type II (Failure) Censoring . . . . . . . . . . . . . . . 38
12.3.2 Type I Censoring . . . . . . . . . . . . . . . . . . . . . 38
13 Degradation Data, Models, and Data Analysis 38
13.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
13.2 Models for Degradation Data . . . . . . . . . . . . . . . . . . 39
13.2.1 Degradation Data . . . . . . . . . . . . . . . . . . . . . 39
13.2.2 Degradation Leading to Failure . . . . . . . . . . . . . 39
13.2.3 Models for Variation in Degradation and Failure Times 41
13.2.4 Limitations of Degradation Data . . . . . . . . . . . . 41
13.2.5 General Degradation Path Model . . . . . . . . . . . . 41
13.2.6 Degradation Model Parameters . . . . . . . . . . . . . 42
13.3 Estimation of Degradation Model Parameters . . . . . . . . . 42
13.4 Models Relating Degradation and Failure . . . . . . . . . . . . 42
13.4.1 Soft Failures: Specified Degradation Level . . . . . . . 42
13.4.2 Hard Failures: Joint Distribution of Degradation and
Failure Level . . . . . . . . . . . . . . . . . . . . . . . 43
13.5 Evaluation of F (t) . . . . . . . . . . . . . . . . . . . . . . . . 43
13.5.1 Analytic Solution for F (t) . . . . . . . . . . . . . . . . 43
13.5.2 Numerical Evaluation of F (t) . . . . . . . . . . . . . . 44
13.5.3 Monte Carlo Evaluation of F (t) . . . . . . . . . . . . . 44
13.6 Estimation of F (t) . . . . . . . . . . . . . . . . . . . . . . . . 44
13.7 Bootstrap Confidence Intervals . . . . . . . . . . . . . . . . . 44
13.8 Comparison with Traditional Failure-Time Analysis . . . . . . 46
13.9 Approximate Degradation Analysis . . . . . . . . . . . . . . . 46
14 Introduction to the Use of Bayesian Methods for Reliability
Data 47
14.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
14.2 Using Bayes Rule to Update Prior Information . . . . . . . . . 47
14.3 Prior Information and Distributions . . . . . . . . . . . . . . . 47
14.3.1 Noninformative(diffuse) Prior Distributions . . . . . . . 47
14.3.2 Using Past Data to Specify a Prior Distribution . . . . 47
14.3.3 Expert Opinion and Eliciting Prior Information . . . . 48
14.4 Numerical Methods for Combining Prior Information with a
Likelihood . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
14.4.1 Simulation-based Methods for Computing the Poste-
rior Distribution of θ . . . . . . . . . . . . . . . . . . . 48
14.4.2 Marginal Posterior Distributions . . . . . . . . . . . . . 48
14.5 Using The Posterior Distribution for Estimation . . . . . . . . 48
14.5.1 Bayesian Point Estimation . . . . . . . . . . . . . . . . 48
4
CONTENTS
14.5.2 Bayesian Interval Estimation . . . . . . . . . . . . . . . 49
14.6 Bayesian Prediction . . . . . . . . . . . . . . . . . . . . . . . . 49
14.6.1 Bayesian Posterior Predictive Distribution . . . . . . . 50
14.6.2 Approximating Posterior Predictive Distribution . . . . 50
14.6.3 Posterior Predictive Distribution for the kth Failure
from a Future Sample of Size m . . . . . . . . . . . . . 50
14.7 Practical Issues in the Application of Bayesian Methods . . . . 50
14.7.1 Comparison Between Bayesian and Likelihood/Frequen-
tist Statistical Methods . . . . . . . . . . . . . . . . . . 50
14.7.2 Caution on the Use of Prior Information . . . . . . . . 51
15 System Reliability Concepts and Methods 51
15.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
15.2 System Structures and System Failure Probability . . . . . . . 52
15.2.1 Time Dependency of System Reliability . . . . . . . . . 52
15.2.2 System with Component in Series . . . . . . . . . . . . 52
15.2.3 System with Components in Parallel . . . . . . . . . . 53
15.2.4 Systems with Components in Series-Parallel . . . . . . 53
15.2.5 Bridge System Structure . . . . . . . . . . . . . . . . . 54
15.2.6 k-Out-of-s System Structures . . . . . . . . . . . . . . 54
15.3 Estimating System Reliability From Component Data . . . . . 54
16 Analysis of Repairable System and Other Recurrence Data 54
16.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
16.1.1 Repairable System Reliability Data and Other Recur-
rence Data . . . . . . . . . . . . . . . . . . . . . . . . . 54
16.1.2 A Nonparametric Model for Recurrence Data . . . . . 55
16.2 Non-Parametric Estimation of The MCF . . . . . . . . . . . . 56
16.2.1 Non-Parametric Model Assumptions . . . . . . . . . . 56
16.2.2 Point Estimate of the MCF . . . . . . . . . . . . . . . 56
16.2.3 Standard Errors and Non-parametric Confidence Inter-
vals for the MCF . . . . . . . . . . . . . . . . . . . . . 57
16.2.4 Adequacy of Normal-Approximation Confidence Inter-
vals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
16.3 Non-Parametric Comparison of Two Samples of Recurrence
Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
16.4 Parametric Models for Recurrence Data . . . . . . . . . . . . . 59
16.4.1 Poisson Process . . . . . . . . . . . . . . . . . . . . . . 59
16.4.2 Homogeneous Poisson Process(HPP) . . . . . . . . . . 59
16.4.3 Non-homogeneous Poisson Process(NHPP) . . . . . . . 59
16.4.4 Renewal Process . . . . . . . . . . . . . . . . . . . . . 60
5
CONTENTS
16.4.5 Superimposed Renewal Process . . . . . . . . . . . . . 60
16.5 Tools for Checking Point-Process Assumptions . . . . . . . . . 61
16.5.1 Tests for Recurrence Rate Trend . . . . . . . . . . . . 61
16.5.2 Test for Independent Interrecurrence Times . . . . . . 62
16.6 Maximum Likelihood Fitting of Poisson Process . . . . . . . . 62
16.6.1 Poisson Process Likelihood . . . . . . . . . . . . . . . . 62
16.6.2 Superimposed Poisson Process Likelihood . . . . . . . 63
16.6.3 ML Estimation for the Power NHPP and Loglinear
NHPP . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
16.6.4 Confidence Intervals for Parameters and Functions of
Parameters . . . . . . . . . . . . . . . . . . . . . . . . 63
16.6.5 Prediction of Future Recurrences with a Poisson Process 63
16.7 Generating Pseudo-Random Realizations from An NHPP Pro-
cess . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
16.7.1 General Approach . . . . . . . . . . . . . . . . . . . . . 63
17 Failure-Time Regression Analysis 64
6
1 RELIABILITY CONCEPTS AND RELIABILITY DATA
1 Reliability Concepts and Reliability Data
1.1 Introduction
Definition of Reliability: the probablity that a system, vehicle, machine
etc will perform its intended function under operating conditions, for a spec-
ified period of time.
Reliability is quality over time.
Features of Reliability Data:
1. Typically censored due to the frequent need to analyze life test data
before all units are failed.
2. Most reliability data are modeled using distributions for positive
random variables such as exponential, Weibull . . .
3. Inferences and predictions involving extrapolation are often required.
4. It’s often necessary to use past experience or other scientific or
engineering judgement to provide infomation as input to the analysis
of data or to a decision makeing process.
5. The traditional parameters (such as mean and standard deviation) of
a statistical model are typically not of primary interest. Instead,
specific measures of product reliability or particular characteristics of
a failure time distribution (eg. quantiles, rates) are of more interest.
6. Model fitting requires computer implementation of numerical
methods, and often there is no exact theory for statistical inferences,
especially with censored data.
1.2 Genral Models for Reliability Data
1.2.1 Target Population or Process
Deming(1975) stated that statistical studies can be broadly divided into
two different categories:
1. Enumerative studies answer questions about populations
that consist of a finite set of identifiable units. Typically, the
study is conducted by randomly selecting a sample from the
7
1 RELIABILITY CONCEPTS AND RELIABILITY DATA
population, carefully evaluating the units in the sample and then
making an inference about the larger population from which the
sample is taken.
Assumption: the sampling frame accurately represents the units in
hte population.
2. Anlytic studies answer questions about processes that
generate units or other output over time. Interest might center
on the life distribution of electric motors that will be produced in the
future.
Assumption: The process will behave in the future the same as it
has in the past.
1.2.2 Causes of Failure and Degradation Leading to Failure
Many failure modes can be traced to some underlying degradation process.
For example, automotive brake pads wear with use.
Traditionaly, most statistical studies of product reliability have been based
on failure time data, while the actual level of degradation can also be very
helpful, especially when there are few or no failures.
Not all failrues can be traced to degradation, though. Some failures can be
caused by sudden accidents.
Understanding the physical and chemical mechanisms and random risk
factors leading to failure can suggest methods for eliminating failure modes
or reducing the probablity of a failure mode, thereby improving reliability.
1.2.3 Environment Effects on Reliability
Environmental factors play an important role in product reliability. A large
proportion of product reliability problems result from unanticipated failure
modes caused by environmental effects that were not part of the initial
reliability evaluation system.
One challenge of product design is to discover and develop economical
means of building in robustness to environmental and other factors that
manufactures and users are unable to control.
8
1 RELIABILITY CONCEPTS AND RELIABILITY DATA
1.2.4 Time Scale
The life of many products can be measured in more than one scale. For
example, life would be measured in number of cycles for factory life tests of
products such as wahsing machines and toasters, while time in service data
is more commonly available.
The choice of a time for measuring product life is often suggested by an
underlying process leading to failure, even if the degradation cannot be
observed directly.
There are a number of methods that can be used to handle data with more
than one life measurements. For example, we measure battery age and
number of charge/discharge cycles. We can 1) estimate the effects of both
factors and develop a measure of battery life as a function of both or 2)
develop a model that use the number of cycles to help explain the
variability in the time to failure.
1.3 Repairable Systems and Nonrepairable Units
There are two situations regarding to failure:
1. The time of failure for nonrepairable units or components, or time to
first failure of a system.
2. A sequence of reported system-failure times for a repairable system.
A system may have both repairable and nonrepairable components.
1.4 Strategy for Data Collection, Modeling and
Analysis
Reliability studies involving laboratory experimentation or field tracking
require carefully planning and excution. Mistakes can be costly in terms of
materials and time, and there is also a possibility of drawing erroneous
conclusion.
1.4.1 Planning a Reliability Study
1. Define problem to be solved and questions to be answered.
2. Consider the resources available for the study(time, money,equipment
etc.).
9
1 RELIABILITY CONCEPTS AND RELIABILITY DATA
3. Design the experiment or study with an assessment of precision of
estimates as a function of the size of the study. As estimation of
precision generally depends on unknow model parameters, assessment
will require planning values of unknown population and process
characteristics.
It’s often useful to conduct a pilot study to obtain information
needed when little is known about the target population of process.
4. In new situations, it’s useful, before the test, to conduct a trial
analysis of data simulated from a proposed model suggested from
available information/previous experience. See more details at
Chapter 10 , 20 and section 22.5.
1.4.2 Strategy for Data Analysis and Modeling
General strategy:
1. Begin the analysis by looking at the data without make any
distributional or other strong model assumptions, i.e., look at the
model without any possible distortion.
The primary tool for this step is graphical analysis.
2. If possible, fit one or more parametric models to the data for the
purpose of description, estimation or prediction.
3. Examine appropriate diagnostics and other use other tools to assess
the adequacy of model assumptions before makeing any estimation or
prediction.
Especially when there is little data, it may be difficult to detect
departures from model assumptions, but this just means that we have
no strong evidence against the assumption but NOT means that the
assumption can be trusted.
4. Proceed with caution to esimate parameters and make predictions if
step 3 is fine. The estimates and predictions usually come with
statistical intervals to reflect uncertainty and variability.
5. Display results of the analysis graphically, including estimates or
predictions with uncertainty bounds.
6. It’s possible to use the results to draw conclusions about the
reliability, contigent on particular model asssumptions.
10
2 MODELS, CENSORING, AND LIKELIHOOD FOR FAILURE-TIME DATA
If assumptions that data provide little information to assess their
adequacy, it’s useful to vary assumptions and assess the impact of
such perturbations have on final answers.
This sensitivity analyses should be reported along with conclusions.
2 Models, Censoring, and Likelihood for
Failure-Time Data
2.1 Models for Continuous Failure-Time Process
The most widely used metric for reliability of a product is its failure-time
distribution and the most commonly collected reliability data.
Cumulative distribution function, probability density function, survival
function and hazard function are used to characterize the probability
distribution for failure time T .
1. Cumulative Density Function(CDF): F (t) = P (T ≤ t) is the cdf
of T. It’s the proportion of the units (will) fail before time t.
2. Probability Density Function(PDF): the derivative of cdf. It
represents the relative frequency of failure times as a function of time.
3. Survival Function(SF): also called as the reliability function, the
complement of the cdf. S(t) = 1 − F (t). It represents the probability
of surviving beyond time t.
4. Hazard Function(HF), also called hazard rate, the instantaneous
failure rate function.
h(t) = lim
∆t→0
P (t < T < t + ∆t|T > t)
∆t
(2.1)
= lim
∆t→0
P (t < T < t + ∆t)
∆t · P (T > t)
(2.2)
=
f(t)
1 − F (t)
(2.3)
[4.1] Cumulative Hazard Function: H(t) =
R
t
0
h(x)d(x), and
F (t) = 1 − exp(−H(t)) = 1 − exp(−
R
t
0
h(x)d(x).
11
2 MODELS, CENSORING, AND LIKELIHOOD FOR FAILURE-TIME DATA
[4.2] Average Hazard Rate: AHR(t
1
, t
2
) =
R
t
2
t
1
h(u)d(u)
t
2
−t
1
=
H(t
2
)−H(t
1
)
t
2
−t
1
,
is the hazard rate over the interval and it’s the approximate fraction
failing per unit time over the specified interval.
[4.3] Hazard rate in FITs: A fit(failure in time) rate is defined as
the hazard function in units of 1/hours, and multiplied by 10
9
, it’s
commonly used in high reliability electronics applications.
12
2 MODELS, CENSORING, AND LIKELIHOOD FOR FAILURE-TIME DATA
Figure 1: The Four Functions
13
2 MODELS, CENSORING, AND LIKELIHOOD FOR FAILURE-TIME DATA
5. Quantile Function:the p quantile of F (t) is the smallest t satisfying
P (T ≤ t) = F (t) ≥ p.
2.2 Models for Discrete Data From a Continuous
Process
1. Partition the time line (0, ∞) into (m+1) observation intervals:
(t
0
, t
1
], (t
1
, t
2
], ··· , (t
m
, t
m+1
], where t
m+1
= ∞andt
0
= 0.
2. Define multinomial probability that a unit will fail in interval i as
π
i
= P (t
i−1
< T ≤ t
i
) = F (t
i
) − F (t
i−1
), where π
i
≥ 0 and
P
m+1
j=1
= 1.
3. Then the conditional probability that a unit will fail in interval i
given that the unit was still operating at the beginning of the interval
i: p
i
= P (t
i−1
< T ≤ t
i
|T > t
i−1
) =
π
i
S(t
i−1
)
.
4. The survival function is then
S(t
i
) =
i
Y
j=1
[1 − p
j
], i = 1, ··· , m + 1. (2.4)
5. The cdf at time t
i
:
F (t
i
) = 1 −
i
Y
j=1
[1 − p
j
] =
i
X
j=1
π
j
. (2.5)
2.3 Censoring
2.3.1 Censoring Mechanism
Type I censoring:remove all unfailed units from the test at a
prespecified time is know as “time censoring” .
Type II censoring: terminate a test after a specified number of failures
is known as “failure censoring”.
Interval censoring: failures are only observed when inspection conducted,
thus interval-censored observations have upper and lower bound on a failure
time, also know as inspection data, grouped data, or read-out data.
Random right censoring: some products may have more than one causes
of failures and if the primary interest in focused on one particular cause,
14
2 MODELS, CENSORING, AND LIKELIHOOD FOR FAILURE-TIME DATA
then the failures from other causes can , in some situations , be viewed as
random right censored data.
Systematic multiple censoring:censoring due to staggered entry of
units. The data may be analyzed when not all units have failure, and the
data will usually be multiply right-censored with some failure times
exceeding some of the running times.
2.3.2 Assumptions on Censoring Mechanism
1. The censoring time of a unit depend ONLY on the history of the
observed failure-time process.
2. Censoring is non-informative, i.e., the censoring times of units provide
no information about the failre-time distribution, this is related to the
first assumption.
2.4 Likelihood
The total likelihood , or joint probability of the data for n independent
observations including left-censored, interval-censored and right-censored
data:
L(p; DAT A) = C ·
m+1
Y
i=1
[F (t
i
)]
`
i
· [F (t
i
) − F (t
i−1
)
d
i
] · [1 − F (t
i
)]
r
i
. (2.6)
where n =
P
m+1
j=1
(d
j
+ `
j
+ r
j
) and C is a constant which is usually taken as
1.
Likelihood for Random Censoring in the Intervals
L
i
(p; data
i
) = {P [(T ≤ C) ∩ (t
i−1
< T ≤ t
i
)]}
d
i
{P [(C ≤ T ) ∩ (t
i−1
< C ≤ t
i
)]}
r
i
=
h
Z
t
i
t
i−1
f
T
(t)[1 − F
C
(t)]dt
i
d
i
×
h
Z
t
i
t
i−1
f
C
(t)[1 − F
T
(t)]dt
i
r
i
,
for r
i
right censored observations, and d
i
failures in (t
i−1
, t
i
].
15
3 NON-PARAMETRIC ESTIMATION
3 Non-parametric Estimation
3.1 Estimation from Singly Censored Interval Data
The non-parametric estimator
ˆ
F (t) based on binomial distribution:
ˆ
F (t) =
#offailuresuptotimet
i
n
=
P
i
j=1
d
j
n
. (3.1)
where n is the initial number of units, d
i
is the units that failed/died in
interval (t
i−1
, t
i
).
3.2 Basic Ideas of Statistical Inference
The estimates of
ˆ
F (t) can be interpreted beyond the particular sample
units, and used to make inference about the process or larger existing
population of units from which the sample units are choosen randomly
(Inferential Statistics).
When repeating the sampling and estimation a lot of times, the sampling
distribution of of
ˆ
F (t
i
) can be produced, which can provide insight into the
possible deviation between
ˆ
F (t
i
) and F (t
i
).
Confidence intervals are one of the most useful ways of quantifying the
uncertainty due to “sampling error” arising from limited sample sizes.
Converage probability is the probability that a confidence interval
procedure will result in an interval containing the quantity of interest.
3.3 Confidence Intervals from Complete or Singly
Censored Data
3.3.1 Point-wise Binomial-based Confidence Interval for F (t
i
)
A conservative 100(1 − α)% confidence interval [F
∼
,
˜
F ] is:
F
∼
(t
i
) =
n
1 +
(n − n
ˆ
F + 1)F
1−α/2;2n−2n
ˆ
F +2,2n
ˆ
F
n
ˆ
F
o
−1
. (3.2)
˜
F (t
i
) =
n
1 +
(n − n
ˆ
F )
(n
ˆ
F + 1)F
1−α/2;2n
ˆ
F +2,2n−2n
ˆ
F
o
−1
. (3.3)
where
ˆ
F is
ˆ
F (t
i
), F is the F distribution.
16
3 NON-PARAMETRIC ESTIMATION
3.3.2 Pointwise Normal-Approximation Confidence Interval for
F (t
i
)
An approximate 100(1 − α)% confidence interval for F (t
i
) is
[F
∼
(t
i
),
˜
F (t
i
)] =
ˆ
F (t
i
) ± z
1−α/2
ˆse
ˆ
F
. (3.4)
where ˆse
ˆ
F
=
q
ˆ
F (t
i
)[1 −
ˆ
F (t
i
)]/n.
3.4 Estimation from Multiply Censored Data
The size of the risk set at the beginning of interval i is:
n
i
= n −
i−1
X
j=0
d
j
−
i−1
X
j=0
r
j
(3.5)
where d
j
is the number of units that died/failed in the j
th
interval (t
i−1
, t
i
],
and r
j
is the number of units censored at t
j
.
Then the estimator of the conditional probability of failing in interval i
given the unit enters this interval:
ˆp
i
=
d
i
n
i
(3.6)
The survival function and non-parametric estimator of F (t
i
) can be
obtained through (2.4).
3.5 Pointwise Confidence Intervals from Multiply
Censored Data
3.5.1 Approximate Variance of
ˆ
F (t
i
)
V ar(
ˆ
F (t
i
)) = V ar(
ˆ
S(t
i
)) as they sum to a constant 1. The approximate
variance can be obtained by first-order Taylor series expansion:
ˆ
S(t
i
) ≈ S(t
i
) +
i
X
j=1
∂S
∂q
j
q
j
(ˆq
j
− q
j
), (3.7)
where q
j
= 1 − p
j
and q
j
values are approximately uncorrelated binomial
proportions.
17
3 NON-PARAMETRIC ESTIMATION
Then
V ar(
ˆ
F (t
i
)) = V ar(
ˆ
S(t
i
)) ≈ [S(t
i
)]
2
i
X
j=1
p
j
n
j
(1 − p
j
)
. (3.8)
Note the detailed derivation of 3.8 is, for a single j
th
term of the
summation, we have
∂S
∂q
j
=
Q
i−1
j=1
q
j
and the variance for this specific term is
(
∂S
∂q
j
)
2
V ar( ˆq
j
) = (S(t
i−1
))
2
q(1−q)
n
= (S(t
i−1
))
2
q
2
(1−q)
nq
= (S(t
i
))
2
1−q
nq
.
3.5.2 Greenwood’s Formula
Substitute p
j
, S(t
i
) by ˆp
j
,
ˆ
S(t
i
) respectively in (3.8) gives the Greenwood’s
formula.
3.5.3 Pointwise Normal-Approximation Confidence Interval for
ˆ
F (t
i
)
A normal approximate 100(1 − α)% confidence interval
[F
∼
,
˜
F ] =
ˆ
F (t
i
) ± z
(1−α/2)
ˆse
ˆ
F
, (3.9)
and it’s based on the assumption that Z
ˆ
F
=
ˆ
F (t
i
)−F (t
i
)
ˆse
ˆ
F
can be approximated
by a standard normal distribution.
When the sample size is not large and the data may be heavily skewed, a
logit transformation may be helpful for a better approximation as
log(
ˆ
F (t
i
)) may behave more like a standard normal random variable, i.e.,
Z
logit(
ˆ
F )
=
logit(
ˆ
F (t
i
)) − logit(F (t
i
))
ˆse
logit(
ˆ
F )
. (3.10)
3.6 Estimation from Multiply Censored Data with
Exact Failures
Exact failures arise from a continuous inspection process, while the width of
intervals approaches zero, the step functions of
ˆ
F (t) increases at the
reported failure times. This limiting case of interval-based non parametric
estimatro is know as the product-limit or Kaplan-Meier estimator.
18
4 LOACTION-SCALE-BASED PARAMETRIC DISTRIBUTIONS
3.7 Simulataneous Confidence Bands
The motivation is to quantify the sampling uncertainty simutaneously over
a range of time t instead of pointwise confidence intervals at particular
specified value of t.
Approximate 100(1 − α)% simultaneous confidence bands for F (t) can be
obtained from
[F
∼
(t),
e
F (t)] =
ˆ
F (t) ± e
(a,b,1−α/2)
ˆse
ˆ
F (t)
, (3.11)
where t ∈ [t
L
(a), t
U
(b)] and the approximate factor e
(a,b,1−α/2)
was computed
from a large-sample approximation given in Nair(1984).
3.8 Uncertain Censoring Times
1. If censoring times are random and the form of the distribution is
know, a likelihood estimation method could be applied.
2. If censoring pattern and distribution unknown, ˆp
i
= d
i
(n
i
− r
i
/2) is
applied.
4 Loaction-Scale-Based Parametric
Distributions
4.1 Quantities of Interest in Reliability Applications
1. Probability of failure p = P (T ≤ t) = F (t; θ) at a specified t.
2. p quantile of the distribution of T with t
p
= F
−1
(p; θ).
3. hazard function with h(t) =
f(t;θ)
1−F (t;θ)
.
4. average/expectation/first moment
E(T ) =
Z
∞
0
tf(t; θ)dt (4.1)
= −
Z
∞
0
td(1 − F (t; θ)) (4.2)
= −tF (t; θ)
∞
0
+
Z
∞
0
(1 − F (t; θ))dt (4.3)
=
Z
∞
0
(1 − F (t; θ))dt (4.4)
19
4 LOACTION-SCALE-BASED PARAMETRIC DISTRIBUTIONS
is a measure of the center of f (t; θ).
One possible confusion is if we rewrite f(t; θ) in terms of the hazard
function above, the equation will be E(T ) =
R
∞
0
[1 − F (t)]th(t)dt.
However, these equations are equivalent: because 1 − F (t) = e
−H(t)
,
H(t) =
R
t
0
h(x)dx and d(e
−H(t)
) = −e
−H(t)
d(H(t)) = −e
−H(t)
h(t)dt,
E(T ) =
Z
∞
0
[1 − F (t)]th(t)dt (4.5)
=
Z
∞
0
e
−H(t)
h(t)tdt (4.6)
= −
Z
∞
0
t d(e
−H(t)
) (4.7)
= −te
−H(t)
∞
0
+
Z
∞
0
e
−H(t)
dt (4.8)
=
Z
∞
0
(1 − F (t)dt (4.9)
It’s the same.
5. variance of T with V ar(T ) =
R
∞
0
[t − ET ]
2
f(t; θ)dt and standard
deviation sd(T ) =
p
V ar(T ), is the spread of the distribution of T .
6. coefficient of variation of T γ
2
= sd(T )/E(T ), is useful for comparing
the relative amount of variability in different distributions. 1γ
2
is the
“signal to noise ratio”.
7. standardized third central moment/ coefficient of skewness with
γ
3
=
R
∞
0
[t−E(T )]
3
f(t;θ)dt
[V ar(T )]
3
2
.
4.2 Location-scale and Log-location-scale
Distributions
Location-scale family distributions can be written as the form:
P (Y ≤ y) = Φ(
y − µ
σ
) (4.10)
While for Log-location-scale family random variable T , log(T) belongs to
the location-scale family.
20
4 LOACTION-SCALE-BASED PARAMETRIC DISTRIBUTIONS
4.3 Exponential Distribution
F (t; θ, γ) = 1 − e
−
t−γ
θ
. (4.11)
Hazard function is constant.
Useful for some kinds of electronic components(for example, capacitors, or
high qualify integrated circuits).
4.4 Normal Distribution
F (y; µ, σ) = Φ
nor
(
y − µ
σ
). (4.12)
Increasing hazard function.
4.5 Log-normal Distribution
F (y; µ, σ) = Φ
nor
(
log(y) − µ
σ
). (4.13)
Appropriate for time to failure caused by a degradation process with
combinations of random rate constants that combine multiplicatively.
Widely used to describe the time to fracture from fatigue crack growth in
metals.
4.6 Smallest Extreme Value Distribution
F (y; µ, σ) = Φ
sev
(
y − µ
σ
). (4.14)
where Φ
sev
= 1 − e
−e
z
in standardized case.
Exponentially increasing hazard function suggests that SEV is suitable for
modeling the life of a product that experiences very rapid wearout after a
certain age.
4.7 Weibull Distribution
F (T ≤ t; η, β) = 1 − e
−(
t
η
)
β
, t > 0 (4.15)
Relationship with SEV: If T has a Weibull distribution, then
Y = log(T ) ∼ SEV (µ, σ) with µ = log(η), σ = 1/β.
21
4 LOACTION-SCALE-BASED PARAMETRIC DISTRIBUTIONS
4.8 Largest Extreme Value Distribution
F (y; µ, σ) = Φ
lev
(
y − µ
σ
). (4.16)
where Φ
lev
= e
−e
−z
in standardized case.
Less commonly used due to the positive probability of negative
observations. It can be used as a model for life if σ is small relative to
µ > 0.
4.9 Logistic Distribution
F (y; µ, σ) = Φ
logis
(
y − µ
σ
). (4.17)
where Φ
logis
=
e
z
1+e
z
in standardized form.
The shape is similar to that of the normal distribution except longer tails,
however, the behavior of the hazard function is different at upper tail.
Preferred to the normal distribution because its cdf can be written in a
simple closed form, if computation is a concern.
4.10 Log-logistic Distribution
F (y; µ, σ) = Φ
logis
(
log(y) − µ
σ
). (4.18)
4.11 Generating Pseudo-random Observations from a
Specified Distribution
4.11.1 Pseudo-random Observations from Continuous
Distributions
Generate U
1
, . . . , U
n
random variables from a UNIF (0, 1), then
t
U
i
= F
−1
t
(U
i
) is a pseudorandom sample from F
T
.
4.11.2 Pseudo-random Observation from Discrete Distributions
Same general idea as Section 4.11.1, but it could be more complicated as
sometimes the discrete quantiles cannot be computed directly.
22
5 OTHER PARAMETRIC DISTRIBUTIONS
4.11.3 Efficient Generation of Censored Pseudorandom Samples
General Approach based on order statistics
U
(i)
denotes the i
th
order statistic from a random sample of size n from
UNIF (0, 1), then
P (U
(i)
≤ u|U
(i−1)
= u
(i−1)
) = 1 −
h
1 − u
1 − u
(i−1)
i
(n−i+1)
, u ≥ u
(i−1)
(4.19)
By applying the method in Section 4.11.1 , generate a random variable U
to replace the P to solve the u as U
(i)
.
Application:
Failure-Censored Samples: generate pseudorandom sample with n units
and r failures.
1. Generate U
1
, . . . , U
(
r
)
(the probability) pseudorandom observations
from UNIF (0, 1).
2. Find the order statistics U
(1)
, . . . , U
(r)
given U
(0)
= 0. by (4.19).
U
(1)
= 1 − [1 − U
(0)
] × (1 − U
1
)
1/n
U
(2)
= 1 − [1 − U
(1)
] × (1 − U
1
)
1/(n−1)
.
.
.
U
(r)
= 1 − [1 − U
(r−1)
] × (1 − U
1
)
1/(n−r+1)
3. The pseudorandom sample from F (t; θ) is T
(i)
= F
−1
[U
(i)
; θ].
Time Censored Sample:
Set t
c
is the cutoff time.
1. Generate a new pseudorandom observation U
i
from UNIF (0, 1).
Compute U
(i)
= 1 − [1 − U
(i−1)
] × (1 − U
i
)
1/(n−i+1)
and
T
(i)
= F
−1
[U
(i)
; θ].
2. If T
i
> t
c
stop; else, i = i + 1 and return to step 1.
5 Other Parametric Distributions
The coefficient of variation (CV) is defined as the ratio of the standard
deviation σ to the mean µ: c
v
=
σ
µ
. It shows the extent of variability in
relation to the mean of the population
1
.
1
https://en.wikipedia.org/wiki/Coefficient_of_variation
23
6 PROBABILITY PLOTTING
6 Probability Plotting
Probability plots use special scales on which a cdf of a particular
distribution plots as a straight line.
6.1 Purpose
1. Assess the adequacy of a particular distributional model.
2. Provide non-parametric graphical estimates of probabilities and
distribution quantiles.
3. Obtain graphical estimation of parametric model parameters by
fitting a straight line through the points on a probability plot.
4. Display the results of a parametric maximum likelihood fit along with
the data.
6.2 Linearing a CDF
By find the transformation of F (t) and t such that the relationship between
the transformed variables is linear, one can obtain a linearized plot of time
versus the CDF.
6.3 Probability Plotting Positions
6.3.1 Criteria for Choosing Plotting Positions
1. Checking distributional assumptions: The choice of plotting positions,
in moderate-to-large samples, is not so important.
2. Estimation of parameters: The “best” plotting positions will depend
on the assumed underlying model and the functions to be estimated.
For complete data, let i index the ordered observations, then p
i
=
i−0.5
n
is a good choice for general purpose use in probability plotting.
3. Display of maximum likelihood fits with data: The criteria is that the
line “fit” the points well when the assumed model being fit with ML
agrees with the data.
24
6 PROBABILITY PLOTTING
6.3.2 Choice of Plotting Positions
1. Continuous inspection data and single censoring:Let t
(i)
be the i
th
ordered failure time,
ˆ
F (t
i
) = i/n be a step function increasing by 1/n
at each reported failure time, a reasonable plotting position is the
midpoint of the jump:
i − 0.5
n
=
1
2
h
ˆ
F (t
(i)
) +
ˆ
F (t
(i−1)
)
i
=
1
2
h
i
n
+
i − 1
n
i
. (6.1)
2. Continuous inspection data and multiple censoring: due to the
multiple censoring, the step increases may be different from 1/n,
p
i
=
1
2
h
ˆ
F (t
(i)
) +
ˆ
F (t
(i−1)
)
i
. (6.2)
3. Interval censored inspection data: Failures are recorded at the upper
endpoint of each interval,
p
i
=
ˆ
F (t
i
). (6.3)
4. Arbitrary censored data: With mixtures of left, right censoring, it
requires a compromise between the other two cases.
6.4 Probability Plots with Specified Shape Parameter
There are some distributions not in the location-scale families and cannot
be transformed into such as distribution. For example, the gamma and
generalized gamma distribution. They may have one or more unknown
shape parameters. Two approaches to specify an unknown shape parameter
for a probability plot:
1. Plot the data with different given values of shape parameter to find a
value that will give probability plot that is nearly linear.
2. Use parametric maximum likelihood methods to estimates the shape
parameter and use the estimated value to construct the probability
plotting scales.
25
7 PARAMETRIC LIKELIHOOD FITTING CONCEPTS: EXPONENTIAL DISTRIBUTION
7 Parametric Likelihood Fitting Concepts:
Exponential Distribution
7.1 Introduction
The appeal of maximum likelihood(ML) stems from the fact that it can be
applied to a wide variety of statistical models and kinds of data(e.g.,
continuous, discrete, categorical, censored, truncated).
Statistical theory shows that, under standard regularity conditions, ML
estimators are “optimal” in large samples(ML estimators are consistent and
asymptotically efficient, i.e., among consistent competitors to ML
estimators, non has a smaller asymptotic variance).
Besides Bayesian methods, there is no general theory that suggests
alternatives to Ml that will be optimal in finite samples.
Statistical modeling, in practice, is an iterative procedure of fitting
proposed models in search of a model that provides an adequate description
of the population or process of interest without being unnecessarily
complicated.
7.2 Parametric Likelihood
Similar with equation (2.6) for non-parametric models, the parametric
likelihood function can be written as a joint probability,
L(θ; DAT A) = C
n
Y
i=1
L
i
(θ; data
i
). (7.1)
In practice, log-likelihood is more widely used with a naturally sums,
L(θ) = log[L(θ)] =
n
X
i=1
L
i
(θ). (7.2)
7.3 Confidence Intervals for θ
7.3.1 Likelihood Confidence Intervals for θ
R(θ) =
L(θ)
L(
ˆ
θ)
is the relative likelihood for θ, then an approximate
100(1 − α)% likelihood-based confidence interval is the set of all values of θ
26
7 PARAMETRIC LIKELIHOOD FITTING CONCEPTS: EXPONENTIAL DISTRIBUTION
such that
− 2log(R(θ)) ≤ χ
2
1−α;n
(7.3)
, where n is the degree of number of parameters in θ.
If −2log(R(θ
0
)) > χ
2
1−α;n
, then the hypothesis testing θ = θ
0
should be
rejected.
7.3.2 Normal-Approximation Confidence Intervals for θ
A 100(1 − α)% normal-approximation confidence interval for θ is
[θ
∼
,
e
θ] =
b
θ ± z
(1−α/2)
ˆse
ˆ
θ
, (7.4)
where ˆse
ˆ
θ
=
r
h
−
d
2
L(θ)
dθ
2
i
−1
θ=
ˆ
θ
.
It’s based on the assumption that Z =
ˆ
θ−θ
ˆse
ˆ
θ
can be approximated by
NOR(0, 1).
An alternative approximate confidence interval for positive quantities is,
[θ
∼
,
e
θ] = [
ˆ
θ/w,
ˆ
θ × w], (7.5)
where w = e
z
(1−α/2)
bse
ˆ
θ
/
ˆ
θ
.
log(
ˆ
θ) ± z
1−α/2
bse
log(
ˆ
θ)
bse
log(
ˆ
θ)
≈ bse
ˆ
θ
/
ˆ
θ
It’s based on the assumption that the distribution
Z
log(
ˆ
θ)
=
log(
ˆ
θ) − log(θ)
bse
log(
ˆ
θ)
(7.6)
can be approximated by a NOR(0, 1) distribution.
7.4 Confidence Intervals for Function of θ
For one parameter distributions, confidence intervals for θ can be translated
directly into confidence intervals for monotone functions of θ.
For models have more than one parameter, a collection of intervals must be
handled differently because the confidence level applies only to the process
of constructing an interval for a single point in time t
e
.
27
7 PARAMETRIC LIKELIHOOD FITTING CONCEPTS: EXPONENTIAL DISTRIBUTION
7.5 Likelihood for Exact Failure Times
7.5.1 Correct Likelihood for Observations Reported as Exact
Failures
Sometimes the failure time are reported as exact times, however, are
actually discrete. Then the “correct likelihood” can be corrected as
L
i
(θ; data
i
) =
Z
ti
t
i−1
f(t; θ)dθ = F (t
i
; θ) − F (t
i−1
; θ). (7.7)
7.5.2 Using Density Approximation for Observations Reported
as Exact Values
For most statistical models, the contribution to the likelihood of the
observations reported as exact values can, for small ∆
i
> 0, be
approximated by
F (t
i
; θ) − F (t
i
− ∆
i
; θ) ≈ f(t
i
; θ)∆
i
. (7.8)
where ∆
i
doesn’t depend on θ. Then the likelihood
L
i
(θ) = f(t
i
; θ) (7.9)
differs (7.8) by a constant scale factor. As long as the approximation in
(7.8) is adequate, the general character of the likelihood is not affected.
(7.9) can be used for the likelihood for an observation with “exact failure
time” at t
i
.
7.6 Data Analysis with No Failures
It’s possible that there will be no failures for data on hight-reliability
components. With zero failures and assumed exponential distribution, a
conservative 100(1 − α)% lower confidence bound on θ is
θ
∼
=
2 · T T T
χ
2
1−α;2
=
T T T
−log(α)
, (7.10)
while χ
2
1−α;2
= −2log(α)
2
.
2
Solve the quantile using the chi-square distribution formula, https://en.wikipedia.
org/wiki/Chi-squared_distribution
28
8 MAXIMUM LIKELIHOOD FOR LOG-LOCATION-SCALE DISTRIBUTIONS
8 Maximum Likelihood for
Log-Location-Scale Distributions
8.1 Likelihood
L(µ, σ) =
n
Y
i=1
[f(y
i
; µ, σ)]
δ
i
[1 − F (y
i
; µ, σ)]
1−δ
i
, (8.1)
where
y =
1 if y
i
is an exact observation
0 right-censored observation.
8.2 Likelihood Confidence Regions and Intervals
8.2.1 Joint Confidence Regions for µ and σ
Using the large-sample χ
2
approximation for the distribution of the
likelihood ratio statistic, an approximate 100(1 − α)% joint confidence
region for two dimensional vector,
R(θ
i
, θ
j
) > exp[−χ
2
1−α;2
/2]. (8.2)
8.2.2 Individual Confidence Intervals for µ and σ
The profile likelihood for µ is
R(µ) = max
σ
L(µ, σ)
L(ˆµ; ˆσ)
, (8.3)
Given any µ value, one can calculated and find the point of highest relative
likelihood by maximizing over σ.
The 100(1 − α)% approximate confidence interval interval is given by
R(µ) > exp[−χ
2
1−α;1
/2]. (8.4)
Confidence interval for σ can be derived similarly.
8.2.3 Likelihood Confidence Intervals for Functions of µ and σ
Due to the invariance property of ML estimators, likelihood based methods
can be applied to make inferences about functions of parameters.
For example, for the p quantile, t
p
= exp[µ + Φ
−1
(p)σ],
R(t
p
) = max
σ
L(t
p
, σ)
L(ˆµ, ˆσ)
,
where we substitute µ by t
p
.
29
9 BOOTSTRAP CONFIDENCE INTERVALS
8.3 Normal Approximation Confidence Intervals
1. First compute the variance-covariance matrix. For location-scale
distribution, the local estimate
c
P
ˆ
θ
of
P
ˆ
θ
is the inverse of the
observed information matrix,
d
X
ˆµ,ˆσ
=
"
d
V ar
ˆµ
d
Cov
ˆµ,ˆσ
d
Cov
ˆµ,ˆσ
d
V ar
ˆσ
#
=
"
−
∂
2
L(µ,σ)
∂µ
2
−
∂
2
L(µ,σ)
∂µ∂σ
−
∂
2
L(µ,σ)
∂µ∂σ
−
∂
2
L(µ,σ)
∂µ
2
#
−1
. (8.5)
The partial second derivative describe the curvature of the log
likelihood at the ML estimate, More curvature in the log likelihood
surface implies a more concentrated likelihood near µ, σ, and this
implies more precision.
2. Then for parameters, bse can be derived for normal confidence interval
approximation.
3. For normal approximation confidence intervals for functions of
parameters, say g = g(µ, σ),one can first apply delta method for the
variance of the function, i.e.,
bse
ˆg
=
q
d
V ar(ˆg) =
"
∂g
∂µ
2
d
V ar(ˆµ) + 2
∂g
∂µ
∂g
∂σ
d
Cov(ˆµ, ˆσ) +
∂g
∂σ
2
d
V ar(ˆσ)
#
,
then find the approximate confidence interval based on it.
4. The intervals can be improved slightly by using t
(
p, ν) instead of z
p
.
8.4 Estimation with Given σ
With such information available, one can provide considerably more
precision from limited data. However, the danger is the given value of σ
may be seriously incorrect, resulting in misleading conclusions.
9 Bootstrap Confidence Intervals
Simulation based intervals provide another important method to obtain
exact or more accurate approximate confidence intervals. It’s expected to
be more accurate than the normal approximation methods and competitive
with the likelihood-based methods.
30
9 BOOTSTRAP CONFIDENCE INTERVALS
9.1 Bootstrap Sampling
General Idea: A confidence interval is judged on the basis of how well the
procedure would perform if it were repeated over and over again. A
confidence interval should not be too wide and the coverage probability
(probability that the interval contains the quantity of interest) should be
equal or close to the nominal coverage probability 1 − α. The idea of
bootstrap sampling is to simulate the repeated sampling process and use
the information from the distribution of appropriate statistics in the
bootstrap samples to compute teh needed confidence interval, reducing the
reliance on large-sample approximations.
When the goal is to compute confidence intervals, the usual
recommendation is to use between B = 2000 to B = 5000 bootstrap
samples. Large B values are recommended for estimating the more extreme
quantiles of the bootstrap distribution that are required for higher
confidence levels.
Sampling Methods
1. Fully “parametric” bootstrap procedure: Simulate each sample of size
n from the assumed parametric distribution,using the ML estimates
computed from the actual data to replace the unknown parameters.
The disadvantage is it requires complete specification of the censoring
process. It’s not a problem with simple censoring such as Type I or
II , but more difficult for complicated systematic or random censoring.
2. “Non-parametric” bootstrap sampling: each sample of size n is
obtained by sampling, with replacement, from the actual data cases
in the original data set.
9.2 Confidence Intervals
Take exponential distribution as an example.
Monte Carlo simulation-based methods can be used to obtain a better
approximation to the distribution of Z
log(
ˆ
θ)
than assuming that
Z
log(
ˆ
θ)
˙∼NOR(0, 1). Here are the steps:
1. Simulate B bootstrap samples of size n;
2. Compute the ML estimates for each bootstrap sample;
31
9 BOOTSTRAP CONFIDENCE INTERVALS
3. For the j-th bootstrap sample, compute Z
log(
ˆ
θ
∗
j
)
=
log(
ˆ
θ
∗
j
)−log(
ˆ
θ)
bse
log(
ˆ
θ
∗
j
)
;
Intervals are based on the distribution of t-like statistics. , so the
method is also called “boostrap-t” method.
4. Intervals for θ based on the assumption that the simulation
distribution of Z
log(
ˆ
θ
∗
)
provides a good approximation of the
distribution of Z
log(
ˆ
θ)
:
[θ
∼
,
e
θ] = [
ˆ
θ/w
∼
,
ˆ
θ × ew], (9.1)
where w
∼
= e
z
log(
ˆ
θ
∗
)
(1−α/2)
bse
ˆ
θ
/
ˆ
θ
, ew = e
z
log(
ˆ
θ
∗
)
(α/2)
bse
ˆ
θ
/
ˆ
θ
and
ˆse
ˆ
θ
=
r
h
−
d
2
L(θ)
dθ
2
i
−1
θ=
ˆ
θ
.
Note that the w here is not symmetric as Section 7.3.2 and z is based
on the quantile of the distribution of Z
log(
ˆ
θ
∗
)
.
5. Intervals for θ based on the assumption that the simulated
distribution of Z
ˆ
θ
∗
provides a good approximation to the distribution
of Z
ˆ
θ
.
[θ
∼
,
e
θ] = [
b
θ − z
ˆ
θ
∗
(1−α/2)
bse
ˆ
θ
,
b
θ − z
ˆ
θ
∗
(α/2)
bse
ˆ
θ
], (9.2)
where the z is the quantile of the distribution of Z
ˆ
θ
∗
.
• Bootstrap confidence intervals for monotone functions of θ can be
obtained by doing the transformation.
• For distributions with more than one parameter, for example, the
Weibull distribution, intervals can be obtained in a similar manner,
i.e., base the bootstrap-t confidence intervals on bootstrap evaluations
of the distributions of the t-like statistics, for example,
Z
ˆµ
, Z
log(ˆσ)
, Z
logit(
ˆ
F (t))
. . .
9.3 Percentile Bootstrap Method
When it’s not easy to compute the standard error for an estimate, the
percentile bootstrap described by Efron and Tibshirani(1993) provides a
simple useful alternative.
1. Suppose B simulations , and there is an estimate of θ,
ˆ
θ
∗
j
for
j = 1, ··· , B.
32
10 PLANNING LIFE TEST
2. The 100(1 − α)% percentile bootstrap interval for θ is
[θ
∼
,
e
θ] = [
b
θ
∗
[l]
,
b
θ
∗
[u]
], (9.3)
where
b
θ
∗
[j]
is the bootstrap sample ordered from smallest to largest,
l = B × (α/2), u = B × (1 − α/2). And l, u would be rounded to the
next lowest and next highest integer, respectively.
The advantage is that the interval doesn’t depend the on the transformation
scale of the parameter(“transformation preerving”). It doens’t work as well
as the previous bootstrap-t method for small samples, though.
10 Planning Life Test
10.1 Introduction
Idea: When some “planning information” about the life distribution is
available, its possible to assess the effect of sample size and test length on
the outcome of a particular test plan. Such planning information is
typically obtained from design specifications, expert opinion, or previous
experience with similar product or materials.
represents a planning value of a population or process quantity.
Simulation of a Proposed Test Plan
1. Use the chosen model and planning values of the distribution to
simulate data from the proposed life test.
2. Analyze the data, perhaps fitting more than one distribution.
3. Assess precision of estimates: This can be done initially by computing
approximate confidence intervals, as done in the real data.
4. Simulate many samples and fit distributions for each, assess the
sample-to-sample differences. Such multiple simulations provide an
assessment of estimation precision.
5. Repeat the simulation-evaluation process with different sample sizes
to gauge the actual sample size and test length requirements to
achieve the desired precision.
6. Repeat the simulation-evaluation process with different input
“planning values” over the range of their uncertainty.
33
10 PLANNING LIFE TEST
In general, simulation is a useful method of assessing variability. To control
the standard deviation of an estimator to a specified degree of precision, it’s
possible to interpolate among simulated values at different sample sizes.
10.2 Approximate Variance of ML Estimators
10.2.1 Basic Large-Sample Approximation
For a model with θ = (θ
i
, . . . , θ
k
), the results hold approximately for large
samples:
1.
ˆ
θ follow a multivariate normal distribution with mean vector θ and
covariance matrix
P
ˆ
θ
.
2. The large sample approximate covariance matrix can be computed
from
P
ˆ
θ
= I
−1
θ
with the Fisher information,
I
θ
= E
h
−
∂
2
L(θ)
∂θ∂θ
0
i
=
n
X
i=1
E
h
−
∂
2
L
i
(θ)
∂θ∂θ
0
i
. (10.1)
In practical problems, interest will center on one or more scalar
functions of the parameters, sya g = g(θ). Then in large samples,
I : Assume the distribution of ˆg = g(
ˆ
θ) can be approximated by a normal
distribution,
ˆg ˙∼NOR(g(θ), Ase(ˆg), (10.2)
where
Ase(ˆg) =
1
√
n
p
V
ˆg
V (ˆg) = nAvar(ˆg)
Avar(ˆg) =
h
∂g(θ)
∂θ
i
0
X
ˆ
θ
h
∂g(θ)
∂θ
i
II : If g(θ) is positive for all θ, then an alternate form is better. Assume the
distribution of log(ˆg) = log(g(
ˆ
θ)) can be approximated by a normal
distribution,
log(ˆg) ˙∼NOR(log(g(θ)), Ase(log(ˆg)), (10.3)
where
Ase(log(ˆg)) =
1
√
n
p
V
log(ˆg)
V (log(ˆg)) = nAvar(log(ˆg))
Avar(log(ˆg)) =
1
g
2
Avar(ˆg)
34
11 PARAMETRIC MAXIMUM LIKELIHOOD: OTHER MODELS
10.3 Sample Size for Unrestricted Functions
Assume normal distribution of ˆg, and the actual confidence
intervalhalf-width is D,
D = z
(1−α/2)
1
√
n
q
b
V
ˆg
(10.4)
n =
z
2
(1−α/2)
V
ˆg
D
2
. (10.5)
Similar way is used for sample sizes for other kinds of distributions.
11 Parametric Maximum Likelihood: Other
Models
11.1 Truncated Data and Truncated Distributions
Censored Data Truncated Data
It occurs when there is a bound
on an observation, that is, lower
bounds for observations censored
on the right, upper bounds for obs
censored on the left, and both up-
per and lower bounds for interval
censored data.
It arises when even the existence
of a potential observation would
be unknown, if its value were to
lie in a certain range.
Table 1: Difference between Censored and Truncated Data
1. Likelihood with Left Truncation:
L
i
(θ) = P(t
L
i
< T
i
≤ t
i
|T
i
> τ
L
i
) =
F (t
i
; θ) − F (t
L
i
; θ)
1 − F (τ
L
i
; θ)
, t
i
> t
L
i
≥ τ
L
i
.
(11.1)
For an observation with exact failure time t
i
,
L
i
(θ) =
f(t
i
; θ)
1 − F (τ
L
i
; θ)
, t
i
> τ
L
i
. (11.2)
2. Likelihood with Right Truncation: similar with the left truncation.
35
12 PREDICTION OF FUTURE RANDOM QUANTILES
3. Likelihood with Right and Left Truncation for an interval censored
observation:
L
i
(θ) = P(t
L
i
< T
i
≤ t
i
|τ
L
≤ T < τ
U
) =
F (t
i
; θ) − F (t
L
i
; θ)
F (τ
U
i
; θ) − F (τ
L
i
; θ)
, τ
L
i
≤ t
L
i
< t
i
≤ τ
U
i
.
(11.3)
11.2 Summary
General guidelines for fitting parametric distributions to skewed data.
1. If data is left-skewed, it’s possible to fit a three-parameter Weibull
distribution and achieve a good fit to the data. In many cases,
however, it will be possible to fit the simpler two-parameter Weibull
or smallest extreme value distributions and get, effectively, the same
results.
2. If the data is approximately symmetric, one can generally fit a
three-parameter Weibull, a three-parameter lognormal model, or
two-parameter versions of the distributions and get a reasonably good
fit. In many cases, however, it will be possible to fit the simpler
two-parameter normal or logistic distributions and get, effectively, the
same results.
3. If the data is right skewed, it’s often to fit a three-parameter Weibull
or lognormal distribution to get a good fit of the data.
4. The use of a threshold parameter can be viewed from two different
directions. Sometimes it might be viewed as a physical parameter —
a time before which probability of failure is zero or a threshold
strength. In other cases, γ is one of several parameters of a curve
being fit to the data.
12 Prediction of Future Random Quantiles
12.1 Probability Prediction Intervals(θ Given)
With a completely specified continuous probability distribution, an exact
100(1 −α)% “probability prediction interval” for a further observation from
F (t; θ) is (ignoring the data),
P I(1 − α) = [T
∼
,
e
T ] = [t
α/2
, t
1−α/2
], (12.1)
36
12 PREDICTION OF FUTURE RANDOM QUANTILES
where t
p
is the p quantile of F (t; θ ). The coverage probability of the
interval in (12.1) is,
P [T ∈ P I(1 − α)] = P [T
∼
≤ T ≤
e
T ] = P [t
α/2
≤ T ≤ t
1−α/2
] = 1 − α (12.2)
for a continuous distribution.
12.2 Statistical Prediction Interval(θ Estimated)
12.2.1 Coverage Probability Concepts
In statistical prediction, the objective is to predict the random quantity T
based on “learning” sample information(denoted by DAT A). Generally
with only sample data, there is uncertainty in the distribution parameters.
There are two kinds of coverage probabilities:
1. For fixed DATA (and thus fixed
ˆ
θ and [T
∼
,
e
T ]), the conditional
coverage probability of a particular interval [T
∼
,
e
T ] is,
CP [P I(1−α)|
ˆ
θ; θ] = P (T
∼
≤ T ≤
e
T |
ˆ
θ; θ) = F (
e
T ; θ)−F (T
∼
; θ). (12.3)
This conditional probability is unknown because F (t; θ) depends on
the unknown θ.
2. From sample to sample, the conditional coverage probability is
random because [T
∼
,
e
T ] depends on
ˆ
θ. The unconditional coverage
probability for the prediction interval procedure is
CP [P I(1 − α)|θ] = P (T
∼
≤ T ≤
e
T |θ) = E
ˆ
θ
{CP [P I(1 − α)|
ˆ
θ, θ]}.
(12.4)
This is the one generally used to describe a prediction interval
procedure. In general, CP [P I(1 − α); θ ] ≈ 1 − α because of the
dependency on the unknown θ.
12.2.2 Naive Method for Computing a Statistical Prediction
Interval
A “naive” prediction interval for continuous T is obtained by substituting
the maximum likelihood estimate for θ into (12.1).
37
13 DEGRADATION DATA, MODELS, AND DATA ANALYSIS
12.3 The (Approximate) Pivotal Method for
Prediction Intervals
12.3.1 Type II (Failure) Censoring
A life test is run until a specified number of failures (r). When T has a
log-location-scale distribution, and the data are complete or Type II
censored, the random variable Z
log(T )
is pivotal,i.e., the distribution
depends only on n and r, not θ.
The prediction interval can be obtained from,
P [ˆµ + z
log(T )
α/2
× ˆσ < log(T ) ≤ ˆµ + z
log(T )
1−α/2
× ˆσ] = 1 − α. (12.5)
Then quantiles z
log(T )
α/2
, z
log(T )
1−α/2
can be obtained from the distribution of
Z
log(T )
, which can be obtained approximately(due only to Monte Carlo
error) by simulating B realizations of Z
log(T
∗
)
=
log(T
∗
)−ˆµ
∗
ˆσ
∗
with the following
steps,
1. Draw a sample of size n from a log-location-scale distribution with
parameters (ˆµ, ˆσ), censored at the r
th
failure.
2. Use the censored sample to compute ML estimates ˆµ
∗
and ˆσ
∗
.
3. Draw an additional single observation T
∗
from the log-location-scale
distribution with parameters (ˆµ, ˆσ).
4. Compute Z
log(T
∗
)
=
log(T
∗
)−ˆµ
∗
ˆσ
∗
.
5. Repeat 1 to 4 B times. Obtain the approximation for the quantiles
z
log(T )
α/2
and z
log(T )
1−α/2
from the empirical distribution of Z
log(T
∗
)
.
12.3.2 Type I Censoring
Z
log(T )
is only approximately pivotal and quantiles of Z
log(T )
depend on
F (t
c
; θ), the unknown expected proportion failing by time t
c
, and the
sample size n. Thus the prediction interval to predict T is approximate.
13 Degradation Data, Models, and Data
Analysis
13.1 Introduction
Design of high-reliability systems generally requires that the individual
system components have extremely high reliability, even after a long
38
13 DEGRADATION DATA, MODELS, AND DATA ANALYSIS
periods of time. With short development times, tests must be conducted
with severe time constrains. Frequently no failures occur during such tests.
A relationship between component failure and amount of degradation
makes it possible to use degradation models and data to make inferences
and predictions about failure time.
13.2 Models for Degradation Data
13.2.1 Degradation Data
Degradation Data:
1. The measurement of physical degradation as a function of time.
2. Actual physical degradation cannot be measured, but measures of
product performance degradation may be available.
Both kinds are generally referred to as “degradation data”. Modeling
performance degradation data may be useful but complicated because
performance may be affected by more than one underlying degradation
process.
Advantages of degradation data:
1. Degradation data can provide considerably more reliability
information than traditional censored failure-time data.
2. Accelerated tests are commonly used obtain reliability tests
information more quickly. Direct observation of the physical
degradation process or some closely related surrogate may allow
direct modeling of the failure causing mechanism, providing more
credible and precise reliability estimates and a firmer basis for
often-needed extrapolation.
13.2.2 Degradation Leading to Failure
Most failures can be traced to an underlying degradation process.
Figure 2 on Page 40 shows different possible shapes of the degradation
curves.
39
13 DEGRADATION DATA, MODELS, AND DATA ANALYSIS
Figure 2: Possible Shapes for Univariate Degradation Curves
40
13 DEGRADATION DATA, MODELS, AND DATA ANALYSIS
13.2.3 Models for Variation in Degradation and Failure Times
There is some degree of variability in all of model factors as well as in
factors that are not in the model. These factors combine to cause
variability in the degradation curves and in failure times.
1. Unit-to-Unit Variability
• Initial Conditions
• Material Properties
• Component Geometry or Dimensions
• Within Unit Variability: often there will be spatial variability in
material properties within a unit(e.g., defects).
2. Variability Due to Operating and Environmental Conditions: For
example, stress and temperature.
13.2.4 Limitations of Degradation Data
1. Physical degradation or performance degradation are natural
properties to measure for many testing process. However, degradation
measurements, often, requires destructive inspection or disruptive
measurement. In such situations, one can obtain only one single
measurement on each unit tested.
2. Also, when degradation data are contaminated with large amounts of
measurement error or when the degradation measure is not closely
related to a failure, the advantages of degradation data is
compromised. For example, when the measurement is on performance
degradation, failures may occur for physical reasons that cannot be
observed directly.
An important but difficult engineering challenge of degradation analysis is
to find variables that are closely related to failure time and develop
methods for accurately measuring these variables.
13.2.5 General Degradation Path Model
The actual degradation data of a particular unit over time is denoted by
D(t), t > 0. The observed sample degradation y
i
j of unit i and time t
j
,
y
ij
= D
ij
+
ij
, i = 1, . . . , n, j = 1, . . . , m
i
, (13.1)
41
13 DEGRADATION DATA, MODELS, AND DATA ANALYSIS
where D
ij
= D(t
ij
, β
1i
, . . . , β
ki
) is the actual path of the unit and
ij
∼ NOR(0, σ
). Typically ,β is a vector of k unknown parameters, same
paths have k = 1, 2, 3, 4 parameters.Some β will be random from unit to
unit.
Assumptions:
1. Random β is independent of the
ij
deviations.
2. σ
is constant. The adequacy of this assumption can be affected by
transforming D(t). If there is an autocorrelation among
ij
, j = 1, . . . , m
i
, one can use a time series model for the residual
term along with appropriate estimation methods.
13.2.6 Degradation Model Parameters
We focus on making inferences about the population or process or
predictions about future units. The underlying model parameters are
θ
β
= (µ
β
, Σ
β
), which is the mean and covariance matrix of random β.
13.3 Estimation of Degradation Model Parameters
The likelihood for the random-parameter degradation model in Section
13.2.5,
L(θ
β
, σ
|DAT A) =
n
Y
i=1
Z
∞
−∞
···
Z
∞
−∞
"
m
i
Y
j=1
1
σ
φ
nor
(ζ
ij
)
#
×
f
β
(β
1i
, . . . , β
ki
; θ
β
)dβ
1i
, . . . , dβ
ki
, (13.2)
where ζ
ij
=
y
ij
−D(t
i
j,β
1i
,...,β
ki
)
σ
and f
β
(β
1i
, . . . , β
ki
; θ
β
) is the multivariate
normal distribution density function.
Parameters needed to be estimated are (µ
β
, Σ
β
, σ
).
13.4 Models Relating Degradation and Failure
13.4.1 Soft Failures: Specified Degradation Level
The failure is defined (in a somewhat arbitrary, but purposeful, manner) at
a specified level of degradation if the products have a gradual loss of
performance. For example, light bulb output.
A fixed D
f
, the critical degradation level for failure, will be used.
42
13 DEGRADATION DATA, MODELS, AND DATA ANALYSIS
13.4.2 Hard Failures: Joint Distribution of Degradation and
Failure Level
The definition of a failure is clear — the product stops working.
With hard failures, failure times will not, in general, correspond exactly
with a particular level of degradation. Instead, the level of degradation at
which failure occurs will be random from unit to unit and even over time.
This could be modeled by using a distribution to describe unit to unit
variability in D
f
or a joint distribution of β and stochastic behavior in D
f
,
where D
f
is the critical degradation level for failure.
13.5 Evaluation of F (t)
A specified model for D
t
and D
f
defines a failure-time distribution. This
distribution can be written as a function of the degradation model
parameters. Suppose that if the degradation level first reaches D
f
at time t,
then the unit fails at time t. Then
P (T ≤ t) = F (t) = F (t; θ
β
) = P (D(t; β
1
, . . . , β
k
) ≥ D
f
). (13.3)
That is to say, for a fixed D
f
, the distribution of T depends on the
distribution of the β
1
, . . . , β
k
, which in turn, depends on the basic path
parameters θ
β
.
For most practical models, especially when D(t) is nonlinear and more than
one β
1
, . . . , β
k
are random, F (t) may not be written in a closed form and
numerical evaluation methods need to be used.
13.5.1 Analytic Solution for F (t)
For some simple path models, F (t) can be expressed as a function of the
basic path parameters of a particular unit,
D(t) = β
1
+ β
2
t (13.4)
, where β
1
is fixed and β
2
∼ LOGNOR(µ, σ), i.e., β
1
is the common initial
amount of degradation of all the test units, and β
2
is the random
degradation rates for different units. Then the F (t),
F (t; β
1
, µ, σ) = P [D(t) ≥ D
f
]
= P [β
2
>
D
f
− β
1
t
]
= 1 − Φ
nor
log(
D
f
−β
1
t
) − µ
σ
, t > 0 (13.5)
43
13 DEGRADATION DATA, MODELS, AND DATA ANALYSIS
13.5.2 Numerical Evaluation of F (t)
When a closed form of F (t) exists, one can evaluate it by direct integration.
13.5.3 Monte Carlo Evaluation of F (t)
Monte Carlo simulation is a versatile method for evaluating F (t). The idea
is: generate a large number of random sample paths from the assumed path
model, using the proportion of path crossing D
f
by time t as an evaluation
of F (t).
Algorithm
1. Generate N simulated realizations of
ˇ
β
1
, . . . ,
ˇ
β
k
of β
1
, . . . , β
k
from a
multivariate normal distribution with mean
ˆ
µ
β
and covariance matrix
ˆ
Σ
β
. N is a large number, for example, 10,000.
2. Compute N simulated failure times corresponding to the N
realizations of
ˇ
β
1
, . . . ,
ˇ
β
k
by solving D(t; β
1i
, . . . , β
ki
) = D
f
.
3. For any desired values of t, compute
F (t) ≈
Number of Simulated F irst Crossing T ime ≤ t
N
, (13.6)
and this is an evaluation of F (t).
The error of this approximation is evaluated by using the binomial
distribution: The standard deviation of the Monte Carlo error in F (t) is
p
F (t)(1 − F (t))/N.
13.6 Estimation of F (t)
1. Closed form exists: The failure time distribution F (t) by substituting
the estimates
ˆ
θ
β
into Equation (13.3).
2. No closed form: Algorithm in Section 13.5.2 and 13.5.3 can be used to
evaluate equation (13.3) at
ˆ
θ
β
.
13.7 Bootstrap Confidence Intervals
Bias-corrected percentile bootstrap method is used for obtain
standard errors for
ˆ
F (t).
1. Use the observed data from the n sample paths to compute the
estimates
ˆ
θ
β
and ˆσ
.
44
13 DEGRADATION DATA, MODELS, AND DATA ANALYSIS
2. Use the algorithm in Section 13.5.2 and 13.5.3(either numeric or
Monte Carlo evaluation) to compute the estimate
ˆ
F (t) at desired
values of t. This the estimated F (t).
3. Generate a large number B(e.g. B = 4000) of bootstrap samples that
mimic the original sample and compute the corresponding bootstrap
estimates
ˆ
F
∗
(t) with the following steps:
• Generate n simulated realizations of the random path
parameters β
∗
1i
, . . . , β
∗
ki
, i = 1, . . . , n from
ˆ
θ
β
.
• Using the same sampling scheme as in the original experiment,
compute n simulated observed paths from Equation (13.1)
y
∗
ij
= D
∗
ij
+
∗
ij
, i = 1, . . . , n, j = 1, . . . , m
i
, (13.7)
up to the planned stopping time t
c
, where
∗
ij
is from the normal
distribution based ˆσ
in Step 1, and D
∗
ij
= D(t
ij
; β
∗
1i
, . . . , β
∗
ki
), t
ij
is the same as the original.
• Use the n simulated paths to estimate parameters of the path
model, giving the bootstrap estimates
ˆ
θ
∗
β
.
• Use the algorithm in Section 13.5.2 and 13.5.3(either numeric or
Monte Carlo evaluation) to compute the estimate
ˆ
F
∗
(t) at
desired values of t.
4. For each desired value of t, the bootstrap confidence interval for F (t)
is calculated using the following steps with bias-corrected
percentile bootstrap:
• Sort the B estimates
ˆ
F
∗
(t)
1
, . . . ,
ˆ
F
∗
(t)
B
in increasing order.
• The lower and upper bounds of point-wise approximate
100(1 − α)% confidence intervals for the distribution function
F (t) according to [Tibshirani(1993)]:
h
F
∼
(t),
e
F (t)
i
=
h
ˆ
F
∗
(t)
(l)
,
ˆ
F
∗
(t)
(u)
i
, (13.8)
where
l = B × Φ
nor
2Φ
−1
nor
(q) + Φ
−1
nor
(α/2)
u = B × Φ
nor
2Φ
−1
nor
(q) + Φ
−1
nor
(1 − α/2)
,
and q is the proportion of the B values of
ˆ
F
∗
(t) that are less
than
ˆ
F (t). q = 0.5 is equivalent to the percentile bootstrap.
45
13 DEGRADATION DATA, MODELS, AND DATA ANALYSIS
13.8 Comparison with Traditional Failure-Time
Analysis
Degradation directly models the relationship degradation and time and
takes account of the censored observation when estimating F (t), so
sometimes it can provide a more reasonable extrapolation on F (t) than the
traditional failure time analysis.
13.9 Approximate Degradation Analysis
Approximate degradation analysis provides an (approximately correct)
alternative method of analyzing degradation data.
Steps:
1. First step: separate analysis for each unit to predict the time at which
the unit will reach the critical degradation level corresponding to
failure.
• For the unit i, use the path model y
ij
= D
ij
+
ij
and the sample
path data (t
i1
, y
i1
), . . . , (t
ik
, y
ik
) to find the conditional ML
estimates
ˆ
β
i
= (
ˆ
β
1i
, . . . ,
ˆ
β
ki
).
• Solve D(t;
ˆ
β
i
) = D
f
for t and call the solution
ˆ
t
i
, i.e., the
estimated failure time for unit i.
2. Second step: use all the pseudo failure times from Step 1 as a
complete sample of failure times to estimate F (t): Repeat the
procedure for each sample path to obtain all n pseudo failure times.
Do a single distribution analysis of data
ˆ
t
1
, . . . ,
ˆ
t
n
.
Potential problems:
• Less appealing when the degradation paths are nonlinear.
• It ignores the prediction error in
ˆ
t and doesn’t account for the
measurement error in the observed sample paths.
• The distribution fitted to the pseudo failure times will not, in general,
correspond to the distribution induced by the degradation model.
• For some applications, the sample paths may no contain enough
information to estimate all parameters at some certain time. Fitting
different models for different sample paths could be necessitated.
Overall, extrapolation into the tails of the failure time distribution may be
more valid with the actual crossing distribution implied by the degradation
model than with the empirically predicted failure times.
46
14 INTRODUCTION TO THE USE OF BAYESIAN METHODS FOR RELIABILITY DATA
14 Introduction to the Use of Bayesian
Methods for Reliability Data
14.1 Introduction
Combination of extensive past experience and physical/chemical theory can
provide prior information to form a framework for inference and decision
making. However, there are, of course, dangers involved in making strong
assumptions about knowledge of model parameters.
14.2 Using Bayes Rule to Update Prior Information
The posterior distribution of a set of parameters,
f(θ|DAT A) =
L(DATA|θ)f(θ)
R
L(DATA|θ)f(θ)dθ
=
R(θ)f (θ)
R
R(θ)f (θ)dθ
, (14.1)
where R(θ) =
L(θ)
L(
ˆ
θ)
is the relative likelihood.
There is usually no closed form of this posterior distribution.
14.3 Prior Information and Distributions
Two main sources of prior information:
1. Expert or other subjective opinion.
2. Past data.
14.3.1 Noninformative(diffuse) Prior Distributions
Non-informative pdfs that are constant over the range of the model
parameters. It’s also called “vague prior” or “diffuse prior”.
14.3.2 Using Past Data to Specify a Prior Distribution
Combining past data with a non-informative prior distribution gives a
posterior pdf that is proportional to the likelihood.
47
14 INTRODUCTION TO THE USE OF BAYESIAN METHODS FOR RELIABILITY DATA
14.3.3 Expert Opinion and Eliciting Prior Information
A general approach is to elicit information about particular quantities(or
parameters) that, from past experience(or data), can be specified
approximately independently. For example, for a high reliability integrated
circuit, a good choice would be a quantile in the lower tail of the
failure-time distribution and the lognormal shape parameter σ.
14.4 Numerical Methods for Combining Prior
Information with a Likelihood
14.4.1 Simulation-based Methods for Computing the Posterior
Distribution of θ
Using a larger number of simulated points provides a better approximation
and the number of points used is limited only by computing equipment and
time constrains.
Algorithm 14.1 Monte Carlo Simulation
1. Generate a random sample , θ
i
, i = 1, . . . , M , from the prior f (θ).
2. Retain the i
th
sample θ
i
, with probability R(θ
i
), the relative
likelihood. Do this by generating a random variable U
i
from uniform
distribution and retain the sample if U
i
≤ R(θ
i
).
14.4.2 Marginal Posterior Distributions
Inferences on individual parameters are obtained by using the marginal
posterior distribution of the parameter of interest.
f(θ
j
|DAT A) =
Z
f(θ|DAT A)dθ
−j
(14.2)
Estimates and confidence intervals for a scalar function of the parameters
g(θ) are obtained by using the marginal posterior pdf f[g(θ|DAT A)] and
cdf F [g(θ|DAT A)], which are approximated by the empirical pdf and cdf of
g(θ
∗
) by simulation, respectively.
14.5 Using The Posterior Distribution for Estimation
14.5.1 Bayesian Point Estimation
Bayesian inference for θ and functions of the parameters g(θ) are entirely
based on, f[θ|DAT A] and f[g(θ|DAT A)] respectively. Given that g(θ) is a
48
14 INTRODUCTION TO THE USE OF BAYESIAN METHODS FOR RELIABILITY DATA
scalar, a common Bayesian estimate of g(θ) is the mean of the posterior
distribution.
ˆg(θ) = E[g(θ)|DAT A] =
Z
g(θ)f(θ|DAT A)dθ. (14.3)
This is an estimate that minimizes the square error loss. Other possible
choices to estimate g(θ) are mode of the posterior pdf(similar to the ML
estimate), median and mean. For example,
ˆg(θ) ≈
1
M
∗
M
∗
X
i=1
g(θ
∗
i
),
is the sample mean.
14.5.2 Bayesian Interval Estimation
A 100(1 − α)% confidence interval(credible interval) for a scalar function
g(θ) is any interval satisfying
Z
eg
g
∼
f[g(θ)|DAT A]dg(θ) = 1 − α. (14.4)
The lower and upper bounds can be chosen in different ways:
1. Combine two 100(1 − α/2)% intervals: it puts equal probability in
each tail, preferred when there is more concern for being incorrect in
one direction than the other.
2. Highest posterior density(HPD): it chooses [g
∼
, bg] which consists all
values of g with f(g|DAT A) > c, where c is a constant such that
(14.4) holds. It’s similar to likelihood-based approximate confidence
intervals, calibrated with a χ
2
quantile.
If there are more than one parameters in θ, a 100(1 − α)% confidence
region will be computed:
CR
B
= {g(θ)|f[g(θ)|DAT A] ≥ c}, (14.5)
where c is chosen such that
R
CR
B
f[g(θ)|DAT A]dg(θ) = 1 − α.
14.6 Bayesian Prediction
Bayesian methods are useful for predicting future event like the failure of a
unit from a specified population or process. It can be predicted by using
the Bayesian posterior predictive distribution.
49
14 INTRODUCTION TO THE USE OF BAYESIAN METHODS FOR RELIABILITY DATA
14.6.1 Bayesian Posterior Predictive Distribution
If X represents a future random variable from f(x|θ), then the posterior
predictive pdf of X:
f(x|DAT A) =
Z
f(x|θ)f(θ|DAT A)dθ = E
θ|DAT A
[f(x|θ)]. (14.6)
14.6.2 Approximating Posterior Predictive Distribution
The Bayesian posterior predictive pdf can be approximated by the average
of the posterior pdf f(x|θ
∗
i
),
f(x|DAT A) ≈
1
M
∗
M
∗
X
i=1
f(x|θ
∗
i
), (14.7)
where θ
∗
i
is sampled using Algorithm 14.1. The CDF can be approximated
similarly by replacing f with F .
14.6.3 Posterior Predictive Distribution for the k th Failure from
a Future Sample of Size m
Assume distribution of time T has a log-location-scale distribution. Let T
(k)
denote the k
th
largest observation for a sample of size m from the
distribution of T . The pdf of T
(k)
conditional on θ,
f(t
(k)
|θ) =
m!
(k − 1)!(m − k)!
×[Φ(ζ)]
k−1
×
1
σt
(k)
φ(ζ) ×[1 −Φ(ζ)]
m−k
, (14.8)
where ζ =
log(t
(k)
)−µ
σ
.
The CDF of T
(k)
conditional on θ is,
P [T
(k)
≤ t
(k)
|θ] = F [t
(k)
|θ] =
m
X
j=k
m!
j!(m − j)!
[Φ(ζ)]
j
× [1 − Φ(ζ)]
m−j
. (14.9)
Combining (14.7), (14.8) and (14.9) provides the pdf and cdf of T
(k)
.
14.7 Practical Issues in the Application of Bayesian
Methods
14.7.1 Comparison Between Bayesian and
Likelihood/Frequentist Statistical Methods
Difference: The manner in which nuisance parameters are handled.
50
15 SYSTEM RELIABILITY CONCEPTS AND METHODS
Bayesian Likelihood
Bayesian interval inference meth-
ods are based on a marginal dis-
tribtion in which nuisance param-
eters have been integrated out.
Nuisance parameters can be max-
imized out, as suggested by large-
sample theory.
Parameter uncertainty can be in-
terpreted in terms of probabilities
from the marginal posterior dis-
tribution.
Confidence intervals based on
likelihood and profile likelihood
functions can be calibrated and
interpreted in terms of repeated
sampling coverage probabilities.
In large samples where the likelihood and posterior are approx-
imately symmetric, Bayesian and likelihood confidence interval
methods give very similar answers when prior information is ap-
proximately uninformative.
Table 2: Comparison of Bayesian and Likelihood Methods
14.7.2 Caution on the Use of Prior Information
Analysts and decision makers must beware of and avoid the use of “wishful
thinking” as prior information. The potential for generating seriously
misleading conclusions is especially high when experimental data will be
limited and the prior distribution will dominate in the final answers.
When using Bayesian statistics, it’s important to do sensitivity analysis
with respect to uncertain inputs to one’s model, i.e., change the prior
distribution assumptions and check the effect that the changes have on final
answers of interest.
15 System Reliability Concepts and
Methods
15.1 Introduction
Assessing and improving system reliability generally requires consideration
of system structure and component reliability. Also, some systems are
replaced upon failure, while many are maintained and/or repaired after
failure. For repairable systems, availability may be the appropriate metric.
This leads to consideration of maintainability and repairability. In general,
availability can be increased by increasing reliability or by improving
51
15 SYSTEM RELIABILITY CONCEPTS AND METHODS
maintainability and repairability.
• Availability: the fraction of time taht a system is available for use.
• Maintainability: improvement of reliability through inspection and/or
preventive maintenance.
• Repairability: characterized by the distribution of time to do a repair.
15.2 System Structures and System Failure
Probability
System failure probability, F
T
(t; θ), is the probability that the system fails
before t. The failure probability of the system is a function of
• time in operation,
• the system structure,
• reliability of system components,
• interconnections, and interfaces.
Complicated system structures can generally be decomposed into
collections of the simpler structures, and the methods for evaluation of
system reliability can be adapted to more complicated structures.
15.2.1 Time Dependency of System Reliability
Time dependency of the survival probability is suppressed in this Chapter.
Then the cdf of a system with s components is
F
T
(θ) = g[F
1
(θ
1
), . . . , F
s
(θ
s
)].
15.2.2 System with Component in Series
For a system with s independent components, the cdf is
F
T
(t) = 1 −
s
Y
i=1
(1 − F
i
)). (15.1)
The system hazard function is the sum of component hazard functions,
h
T
(t) =
s
X
i=1
h
i
(t), (15.2)
52
15 SYSTEM RELIABILITY CONCEPTS AND METHODS
which can be derived by hazard function definition in Chapter 2:
F (t) = 1 − e
−H(t)
=
s
Y
i=1
(1 − F
i
) ⇒
H(t) = −log(1 − F (t)) = −
s
X
i=1
log(1 − F
i
) =
X
H
i
.
Importance of Part Count in Product Design:
Rule of thumb in reliability engineering:keep the part count small- keep the
number of individual components in a system to a minimum.
Effect of Positive Dependency in a Two-Component System:
The assumption of independence is conservative in the sense that the actual
F
T
(t) is smaller than that predicted by the independent component model.
If the correlation is negative, the prediction will be anti-conservative. But
negative correlation is uncommon in physical systems.
15.2.3 System with Components in Parallel
For a system with s independent components, the cdf of the system is :
F
T
(t) =
s
Y
i=1
F
i
. (15.3)
Effect of Positive Dependency in a Two-Component
Parallel-Redundant System:
The advantage of redundancy can be degraded seriously when the failure
times of the individual components have positive dependence, as the system
reliability will be smaller if there exits a positive relationship.
15.2.4 Systems with Components in Series-Parallel
Series-Parallel System Structure with System-Level Redundancy
A r × k series-parallel system-level redundancy structure has r parallel
setes, each with k components in series. With independent components, the
cdf is
F
T
(t) =
r
Y
i=1
1 −
k
Y
j=1
(1 − F
ij
)
. (15.4)
Series-Parallel System with Component-Level Redundancy A k ×r
series-parallel system with independent components has k series structures
53
16 ANALYSIS OF REPAIRABLE SYSTEM AND OTHER RECURRENCE DATA
and r components inin parallel.It has cdf ,
F
t
(t) = 1 −
k
Y
i=1
1 −
r
Y
j=1
F
j
. (15.5)
15.2.5 Bridge System Structure
Depending on if the component at the bridge position works or not, the
bridge system can be treated as series-parallel with component level
redundancy or system level redundancy, so the cdf can be computed by the
sum of two conditional cases.
15.2.6 k-Out-of-s System Structures
A system works if at least k out of s components work but not otherwise.
The cdf ,
F
t
(t) = P (T ≤ t)
= P (atleast(s − k + 1)componentsf ail)
=
s
X
j=s−k+1
n
X
δ∈A
j
s
Y
i=1
F
δ
i
i
(1 − F
1−δ
i
i
)
o
.
15.3 Estimating System Reliability From Component
Data
Maximum likelihood estimation for each component and then the mle of
the system can be obtained accordingly. Delta method can be used to find
the variance-covariance matrix. Normal approximation confidence intervals
or bootstrap approximate confidence intervals are proper.
16 Analysis of Repairable System and
Other Recurrence Data
16.1 Introduction
16.1.1 Repairable System Reliability Data and Other
Recurrence Data
1. Generally, repair times are measured in terms of a system age or time
since some well-defined specific event in the system’s history.
54
16 ANALYSIS OF REPAIRABLE SYSTEM AND OTHER RECURRENCE DATA
2. The stochastic model for recurrence data is called a “point-process”
model.
Repairable system data are collected to estimate or predict quantities like:
1. The distribution of times between repairs, τ
j
= T
j
− T
j−1
.
2. The cumulative number of repairs in the interval (0, t] as a function of
system age t.
3. The expected time between failures(MTBF).
4. The expected number of repairs in the interval (0, t] as a function of t.
5. The repair rate as a function of t.
6. Average repair cost as a function of t.
16.1.2 A Nonparametric Model for Recurrence Data
For a single system, recurrence data can be expressed as N(s, t), in system
age interval (s, t], where N(s, t) is the cumulative number of
recurrences in the system. N(t) is used to represent N(0, t) for
simplicity.
1. The model to describe a population of systems is based on the mean
cumulative function(MCF) at system age t.
2. The population MCF is defined as µ(t) = E[N(t)], where the
expectation is over the variability of each system and the unit-to-unit
variability in the population.
3. ν(t) =
dE[N (t)]
dt
=
dµ
dt
, is the recurrence rate per system for the
population given that µ(t) is differentiable.
4. The method introduced next can also be used to model other
quantities accumulating in time than number of repairs, such as the
mean cumulative cost per system.
55
16 ANALYSIS OF REPAIRABLE SYSTEM AND OTHER RECURRENCE DATA
16.2 Non-Parametric Estimation of The MCF
16.2.1 Non-Parametric Model Assumptions
If an observed collection of n ≥ 1 systems is an entire population of interest
or a sample from a larger population of systems, then the method described
here can be used to estimate the population MCF.
1. There exits a population of cumulative functions(one for each system
in the population), from which a sample has been observed.
2. Randomness in the sample is due to the random sampling of
cumulative functions from the population.
3. The time that observation of a system is terminated doesn’t depend
on the system’s history.
Biased MCF may be caused by:
1. Units follow a staggered scheme of entry into service and the
recurrence rate ν(t) is increasing in real time due to some possible
external events affecting all systems simultaneously. Result: The
newer systems that have a more stressful life will be censored earlier
and an overly optimistic estimate of the recurrence rate will be
obtained. The third assumption is not satisfied.
2. Non-parametric estimator doesn’t require that the sampled systems
be statistically independent.
16.2.2 Point Estimate of the MCF
A simple naive estimator of the population MCF:
The sample mean of the available N
i
(t) for the systems still operating at
time t.
• N
i
(t) denotes the cumulative number of system recurrence for system
i before time t.
• Limitation: this estimator is only appropriate if all systems are still
operating at time t.
An appropriate unbiased estimator proposed by Nelson(1988), allows for
different lengths of observation among systems.
56
16 ANALYSIS OF REPAIRABLE SYSTEM AND OTHER RECURRENCE DATA
Algorithm 16.1: Computation of the MCF estimate
1. Order the unique recurrent times, t
ij
among all n systems, where t
ij
denotes the j
th
recurrent time for system i. Let m denote the number
of unique times: t
1
< t
2
< ··· < t
m
.
2. Compute d
i
(t
k
), the total number of recurrences for system i at t
k
.
3. Let δ
i
(t
k
) = 1 if system i is still being observed at time t
k
, and 0
otherwise.
4. Compute
ˆµ(t
j
) =
j
X
k=1
h
P
n
i=1
δ
i
(t
k
)d
i
(t
k
)
P
n
i=1
δ
i
(t
k
)
i
=
j
X
k=1
d
.
(t
k
)
δ
.
(t
k
)
, (16.1)
where δ
.
(t
k
) is the total number of system recurrences at time t
k
and
δ
.
(t
k
) is the size of the risk set at t
k
.
One may plot ˆµ(t) as a piecewise linear function for better visual
perceptions of shape.
16.2.3 Standard Errors and Non-parametric Confidence
Intervals for the MCF
In Nelson(1995a), the true variance of µ(t
j
) for a large population of
cumulative functions is:
V ar[ˆµ(t
j
)] =
j
X
k=1
V ar[
¯
d(t
k
)] + 2
j−1
X
k=1
j
X
ν=k+1
Cov[
¯
d(t
k
),
¯
d(t
ν
)]
=
j
X
k=1
V ar[d(t
k
)]
δ
.
(t
k
)
+ 2
j−1
X
k=1
j
X
ν=k+1
Cov[d(t
k
), d(t
ν
)])
δ
.
(t
k
)
, (16.2)
where
¯
d(t
k
) =
d
.
(t
k
)
δ
.
(t
k
)
is the average number of recurrences per system at t
k
.
Assume that d
i
(t
k
) is a random sample from the population of d(t
k
), and
use moment estimation to obtain estimators of the variance and covariance
above. More details on Chapter 16.2.3 of Meeker and Escobar(1998).
Pointwise normal-approximation confidence intervals for population MCF
can be computed using the general approach and estimated standard error.
57
16 ANALYSIS OF REPAIRABLE SYSTEM AND OTHER RECURRENCE DATA
When µ(t) is positive, intervals can be generated based on
[log(ˆµ(t)) − log(µ(t))]/ ˆse
ˆµ(t)
∼ N(0, 1).
Finite Population Correction: When the number of cumulative
functions sampled is more than 5% or 10% of the population, finite
population methods should be used for estimating standard errors. In this
case,the moment based estimation of variances and covariances should be
multiplied by a factor 1 −
δ
.
(t
k
)
N
, where N is the total number of cumulative
functions in the population of interest.
16.2.4 Adequacy of Normal-Approximation Confidence Intervals
The adequacy depends on:
1. The number of sample cumulative functions at risk to failure;
2. The shape of the distribution of the cumulative function levels at the
point in time where the interval is to be constructed.
Improvements:
1. Use t
p;ν
instead of z
p
, especially when the number of sample systems
at risk is small(say, less than 30).
2. If the cumulative function at a point in time has a normal distribution
and if all units are still under observation, then using t
p;ν
and
[n/(n − 1)]
1/2
se
ˆµ
provide an exact interval for two or more systems.
16.3 Non-Parametric Comparison of Two Samples of
Recurrence Data
Suppose there are two populations/processes with mean cumulative
function µ
1
(t) and µ
2
(t), respectively. Let ∆
µ
(t) = µ
1
(t) − µ
2
(t) represents
the difference at time t.
1. A non-parametric estimator of ∆
µ
(t) is
b
∆
µ
(t) = ˆµ
1
(t) − ˆµ
2
(t).
2. If ˆµ
1
(t) and ˆµ
2
(t) are independent, an estimate of V ar(
b
∆
µ
(t)) is
d
V ar(
b
∆
µ
(t)) =
d
V ar(bµ
1
(t)) +
d
V ar(bµ
2
(t)) .
3. Approximate 100(1 − α)% confidence interval is based on
Z
b
∆
µ
∼ N(0, 1).
58
16 ANALYSIS OF REPAIRABLE SYSTEM AND OTHER RECURRENCE DATA
16.4 Parametric Models for Recurrence Data
The most commonly used models for recurrence data are:
1. Poisson process (homogeneous and non-homogeneous).
2. Renewal process.
3. Superimposed versions of above processes.
16.4.1 Poisson Process
An integer-valued point process on [0, ∞) is said to be a Poisson process
if it satisfies:
1. N(0) = 0.
2. The numbers of recurrences in disjoint time intervals are statistically
independent.
3. The process recurrence rate ν(t) is positive and
µ(a, b) = E[N(a, b)] =
R
b
a
ν(u)du < ∞, when 0 ≤ a < b < ∞.
It follows that N(a, b) has a Poisson distribution:
P [N(a, b) = d] =
µ(a, b)
d
d!
e
−µ(a,b)
d = 0, 1, 2, . . . (16.3)
16.4.2 Homogeneous Poisson Process(HPP)
HPP is a Poisson process with a constant recurrence rate, ν(t) =
1
θ
.
1. The inter-recurrence times, τ
j
= T
j
− T
j−1
, are iid, each with an
exp(θ) distribution, which is followed by
P (τ
j
> t) = p[N(T
j−1
, T
j−1
+ t) = 0].
2. The time T
k
= τ
1
+ ···+ τ
k
to the k
th
recurrence has a Gam(θ, k)
distribution.
16.4.3 Non-homogeneous Poisson Process(NHPP)
NHPP has a nonconstant recurrence rate ν(t). In this case, the
interrecurrence times are neither independent nor identically distributed.
The expected number of recurrences per unit time over (a, b] is:
µ(a, b)
b − a
=
1
b − a
Z
b
a
ν(t)dt. (16.4)
59
16 ANALYSIS OF REPAIRABLE SYSTEM AND OTHER RECURRENCE DATA
An NHPP model is often specified in terms of recurrence rate,
ν(t) = ν(t; θ).
Two examples:
1. Power model recurrence rate, v(t; β, η) =
β
η
(
t
η
)
β−1
, β > 0, η > 0 with
µ(t; β, η) = (
t
η
)
β
.
2. Log-linear model recurrence rate, ν(t; γ
0
, γ
1
) = exp(γ
0
+ γ
1
t) with
µ(t; γ
0
, γ
1
) = [(exp)(γ
0
)][exp(γ
1
t) − 1]/γ
1
.
16.4.4 Renewal Process
A sequence of system recurrences at system ages T
1
, T
2
, ··· is a renewal
process if the iterrecurrence times τ
j
are iid.
The MCF for a renewal process is known as “renewal function”. HPP is a
renewal process.
Before using a renewal process, one needs to check for departures from the
model:
1. Trend and non-independence of interrecurrence times.
2. If the above assumptions are satisfied, one may use earlier chapters’
methods to describe the distribution of interrecurrence times.
16.4.5 Superimposed Renewal Process
The aggregation of renewals from a group of n independent renewal
processes operating simultaneously is known as a superimposed renewal
process(SRP).
1. Unless the individual renewal processes are HPP, an SRP is NOT a
renewal process.
2. When the number of systems n is large and the systems have run long
enough to eliminate transients, an SRP behaves as an HPP. (Similar
with the central limit theorem).
3. The above result on large number of systems can sometimes be used
to justify the use of the exponential distribution to model
interrecurrence times.
60
16 ANALYSIS OF REPAIRABLE SYSTEM AND OTHER RECURRENCE DATA
16.5 Tools for Checking Point-Process Assumptions
16.5.1 Tests for Recurrence Rate Trend
Plots:
1. For a single system, plot of cumulative number of system recurrences
vs. time is the simplest plot. Non-linearity in this plot indicates that
the interrecurrence times are not identically distributed.
2. A plot of interrecurrence times vs. system age, or a “time series plot”
allow the discovery of trends or cycles that would suggest that the
interrecurrence times are not identically distributed.
Tests:
1. Power model: “Military Handbook” test for β = 1(HPP):
X
2
MHB
= −a
r
X
j=1
log(t
j
/t
end
), (16.5)
where t
j
is the t
th
recurrence time, and t
end
is the end of the
observation.
The statistic is approximately χ
2
(2r)
under HPP.
See http://www.itl.nist.gov/div898/handbook/apr/section2/
apr234.htm#TheMilitaryHandbook for more details.
2. Log linear model: Laplace test for trend in log-linear NHPP model,
i.e., γ
1
= 0,
Z
LP
=
P
r
j=1
t
j
/t
end
− r/2
p
r/12
, (16.6)
Z
LP
has approximate standard normally distribution.
3. Both above tests can be misleading when there is no trend but the
underlying process is a renewal process other than HPP. The
Lewis-Robinson test for trend,
Z
LR
= Z
LP
×
¯τ
S
τ
, (16.7)
where ¯τ is the sample mean and S
tau
is the standard deviation of the
inter-recurrence times. The second term on the right hand side is the
reciprocal of the sample coefficient of variation.
In large samples, Z
LR
follows approximately standard normal
distribution under a renewal process. Z
LR
is preferred as a general
test of trend in point-process data.
61
16 ANALYSIS OF REPAIRABLE SYSTEM AND OTHER RECURRENCE DATA
16.5.2 Test for Independent Interrecurrence Times
Consider the serial correlation in the sequence of inter-recurrence times.
Plots: inter-recurrence times vs. lagged inter-recurrence times provides a
graphical representation of serial correlation.
Serial correlation coefficient:
ρ
k
=
Cov(τ
j
, τ
j+k
)
V ar(τ
j
)
. (16.8)
Sample serial correlation:
ˆρ
k
=
P
r−k
j=1
(τ
j
− ¯τ)(τ
j+k
− ¯τ)
q
P
r−k
j=1
(τ
j
− ¯τ)
2
P
r−k
j=1
(τ
j+k
− ¯τ)
2
, (16.9)
where ¯τ =
P
r
j=1
τ
j
/r.
√
r − k × ˆρ
k
˙∼NOR(0, 1) when ρ
k
= 0 and r is large, and this approximate
distribution can be used to asses if ρ
k
is different from zero.
16.6 Maximum Likelihood Fitting of Poisson Process
16.6.1 Poisson Process Likelihood
For ONE system in period (0, t
a
], we have the data as the number of
recurrences, d
1
, ··· , d
m
in non-overlapping intervals (t
0
, t
1
], ··· , (t
m−1
, t
m
]
with t
0
= 0 and t
m
= t
a
. The likelihood for the NHPP according (16.3) is,
L(θ) = P [N(t
0
, t
1
) = d
1
, ··· , N(t
m−1
, t
m
) = d
m
]
=
m
Y
j=1
P [N(t
j−1
, t
j
) = d
j
]
=
m
Y
j=1
[µ(t
j−1
, t
j
; θ)]
d
j
d
j
!
exp
−µ(t
j−1
,t
j
;θ)
=
m
Y
j=1
[µ(t
j−1
, t
j
; θ)]
d
j
d
j
!
× exp
−µ(t
0
,t
a
;θ)
. (16.10)
When the exact reported recurrence times at t
1
≤ ··· ≤ t
r
(r =
P
m
j=1
d
j
),
the likelihood in terms of the density approximation is ,
L(θ) =
m
Y
j=1
ν(t
j
; θ) × exp[−µ(0, t
a
; θ)]. (16.11)
62
16 ANALYSIS OF REPAIRABLE SYSTEM AND OTHER RECURRENCE DATA
16.6.2 Superimposed Poisson Process Likelihood
Assumption: all systems have the same ν(t) (a strong assumption). And
systems are independent.
L(θ) =
n
Y
i=1
[
r
Y
j=1
ν(t
ij
; θ) × exp[−µ(0, t
a
i
; θ)]]. (16.12)
As the assumption is strong, it’s often inappropriate in practice, and
generalizations such as use of explanatory variables to account for
system-to-system differences, are possible.
16.6.3 ML Estimation for the Power NHPP and Loglinear
NHPP
Plug in the ν(t; θ ) into the likelihood function and solve for the mle and
relative likelihood.
16.6.4 Confidence Intervals for Parameters and Functions of
Parameters
General ideas in Chapter 7 and 8 can be used.
16.6.5 Prediction of Future Recurrences with a Poisson Process
The expected number of recurrences in an interval [a, b) is
R
b
a
ν(u; θ)du and
predictions can be made by replacing θ by the maximum likelihood
estimator
ˆ
θ.
16.7 Generating Pseudo-Random Realizations from
An NHPP Process
Simulation data can be used to check the adequacy of large sample
approximation and for implementing bootstrap methods.
16.7.1 General Approach
For a monotone increasing µ(t), the random variable µ(T
i−1
, T
i
), i = 1, ···,
are iid, each with exp(1). The idea is to generate random inter-recurrence
times and solve for random T values sequentially.
1. Generate U
i
, i = 1, ··· , r from a Uniform(0, 1)
63
17 FAILURE-TIME REGRESSION ANALYSIS
µ(T
1
) = −log(U
1
)
µ(T
2
) − µ(T
1
) = −log(U
2
)
.
.
.
µ(T
r
) − µ(T
r−1
) = −log(U
r
). (16.13)
If a realization in an interval (0, t
a
] is needed, then r is random and the
sequential process is stopped when T
i
> t
a
.
General solution:
T
j
= µ
−1
−
j
X
i=1
log(U
i
)
. (16.14)
Or alternatively,
T
j
= µ
−1
µ(T
j−1
) − log(U
j
)
, (16.15)
where T
0
= 0.
17 Failure-Time Regression Analysis
64
