Introduction to Reliability Analysis with Mathematica

Dave Collins
dcollins@outbacksoftware.com

Copyright Outback Software, Ltd., 2002. This notebook may be freely copied and distributed, provided it is kept in its original form, including this notice. For other rights, contact the author. Mathematica is a registered trademark of Wolfram Research, Inc.

Acknowledgements
  
This notebook started as a project for a Mathematica course taught by Ken Levasseur at the University of Massachusetts. Thanks also to Stan Wagon for many hints and tips given at Rocky Mountain Mathematica.

Table of Contents

Introduction

    Initialization

    Mathematical models

    Relibility modeling

Empirical models of reliability

    Non-Parametric estimation of f(t), F(t), R(t) and h(t)

    Example: The exponential distribution

    Maximum likelihood parameter estimation

    Example: The Weibull distribution

    Testing goodness of fit

Markov models

    Reliability modeling examples

Simulation or Monte Carlo methods

    Generating values of random variables

    Reliability modeling example

Summary

Bibliography

Introduction

Reliability analysis is an engineering discipline that applies various mathematical techniques to the measurement and prediction of the reliability of components and systems. The components under study may be mechanical, electronic, software, or other types. Systems could include anything from computers to rail transit. Measurements include failure rates, cumulative failures, and component lifetimes (time until failure). A variety of techniques are employed, drawn mainly from probability, statistics, and the theory of stochastic processes.

Mathematica is a software application for numeric and symbolic computation, developed by Wolfram Research. This being an introduction, the assumption throughout is that the reader knows something about Mathematica, but is not an expert, either in Mathematica or reliability analysis.

The mathematics of reliability analysis, even in terms of specific models, has broad applicability. A reliability model with failures and repairs to failed units may have the same form as a population model in which "failures" become deaths, and "repairs" become births. A probability distribution of survival times may be used in reliability analysis to model component lifetimes, in medical research to model patient survival in a treatment group, and by actuaries to compute insurance premiums.

In this paper a context is established by discussing the general idea of using mathematics to model real-world processes, and outlining the specific methods used by reliability engineers. Following the introduction, several techniques of reliability analysis are discussed in more detail, with particular attention to the suitability of Mathematica as a tool for reliability analysis. Examples have been chosen that use mathematical models with broad applicability, and that provide good illustrations of the use of Mathematica.   

This section consolidates data input and definition of common functions used in this notebook. Most of the cells in this section are Mathematica "initialization cells," meaning they will be automatically evaluated whenever any computational cell in the notebook is evaluated. If you are presented with a dialogue box asking if you want to perform this automatic initialization, always answer "Yes."

The following cell initializes data sets used in subsequent sections. observedPopulation is population data (in millions) for the United States from 1770-2000; ballBearingFailures is time to failure (measured in millions of revolutions) for a set of tested ball bearings; leukemiaRemission is time to spontaneous remission of leukemia in a group of control subjects for a drug trial.   

More typically, data sets like the above would be read in from external files. For example, suppose there is a text file, USpop.txt, wherein each line is a data pair like "1790    3.9", representing the U. S. population (in millions) for a given year. Then the variable observedPopulation could be created as follows:
  
f=OpenRead["USpop.txt"];
observedPopulation=ReadList["USpop.txt",{Number,Number}];
Close[f];
  
This code assumes that USpop.txt is in the default Mathematica folder, as given by Directory[ ]. Since the location of the default directory may vary depending on operating system type and installation options, this is error-prone. Thus the data is imbedded directly in the notebook, for portability. (This would be impractical, of course, for large data sets.)

Data can also be imported in popular spreadsheet formats using built-in Mathematica capabilities. Using MathLink, Mathematica's method of linking to programs written in other language such as C or Java, data could also be read in from external databases such as DB2 or Oracle.
  
Load needed packages:

Turn off unnecessary warning messages:

stepPlot[ ] plots a list of points as a step function (like Histogram[ ], but only showing the outline of the top of the bars). It is useful for plotting things like discrete probability distributions.

It takes one option, PlotWith->{, . . .}, where the fs are other functions to be plotted. The other functions are plotted in a different color. This is useful in comparing the step plot to a fitted function, for example:

(23/(1+ 50Exp[-0.3(#)])&  is the Mathematica "pure function" corresponding to .)

Mathematical models

A mathematical model of a process in the real world is a set of equations with independent variables taking values based on quantitative aspects of the process state at some time, say  , and dependent variables yielding quantitative information regarding the process state at some later time t:

    S(t) = ((t), (t), . . ., (t)) =
        f (S(), , t) = (((), (), . . ., ()), . . ., ((), (), . . ., ()), , t)  

More generally, particularly if the modeling is done on a computer, the function f  may be implemented as an algorithmic procedure that is not necessarily expressible in closed form. A mathematical model may be descriptive (describing or explaining observed data), predictive (predicting data to be observed in the future), or both.

As an example, suppose we are interested in explaining or predicting population growth. One model of growth is that, in an environment with no constraints, future population may be predicted based on a knowledge of the initial population size, and the net growth rate (birth rate minus death rate). It is well known that this model implies exponential growth, as is shown in the function ExponentialGrowthModel below. P0 is the initial population size, and g is the growth rate. In reality, there are always limits to growth, and the logistic model takes this into account, adding a parameter C, the "carrying capacity" of a given environment, which is the maximum size that a population can achieve in that environment before resource limitations reduce the growth rate to zero. The result is a model that first grows exponentially, but then levels off and approaches the carrying capacity asymptotically. This is shown in the function LogisticGrowthModel below.

These two functions are used below to explain population growth in the United States over a two-hundred year period. The actual population figures, at ten-year intervals from census data, are in the ObservedPopulation list (in millions, plotted in black). The parameters used in the two models were estimated visually to produce a reasonable fit (more rigorous techniques of estimating and validating model parameters are described below in Empirical models of reliability). The plot shows that LogisticGrowthModel (plotted in blue) explains the observed data better than ExponentialGrowthModel (plotted in red). (More rigorous methods of measuring "better" are also described below.)

These models are described in [23] and [8]. The 1790-1960 census data is from [9], the 1970-1990 data from [8], and the 2000 data from [24].

There are obvious deficiencies in this exercise. Both models assume a constant growth rate, whereas the actual growth rate over the period is known to have varied. The logistic model, which is a good fit, assumes a constant carrying capacity, and suggestes that the limiting population of the U. S. is about 350 million. In models of animal populations in a fixed territory, availability of food is the main limiting factor, and clearly this has changed over the history of the United States. Possibly there are other factors, such as crowding, that will indeed limit the population as the model suggests, but this remains to be shown. (Failure to recognize the elasticity of human resource limitations was the fallacy of the economist Thomas Malthus, who predicted the "crash" of human population at a level far below what has actually been achieved.)

The choice of a mathematical model involves a tension between several factors:

The model must capture essential aspects of the process under study.

There must be accurate, repeatable measurement procedures for the input and output variables, so that the model can be effectively compared to the real process.

The model must be tractable, both mathematically and economically--i.e., the cost of modelling to a given level of precision must be justified by the value of the result.

The standard advice for building a tractable model is summarized by Friedman ([6], p. 27): "Faced with a complicated problem, assume away any feature that is not essential to what you are trying to understand. When you are finished, you are left with the simplest problem whose solution will tell you what you want to know." In a certain sense this advice is circular--lacking a complete understanding of a problem, we generally can't determine what is essential and what is not. In practice, modeling proceeds iteratively through informed guesswork to trial models which are tested against empirical reality and used to produce more refined models. In spite of the apparently haphazard nature of this process, in the hands of creative scientists it has produced remarkable results. The physicist Eugene Wigner famously remarked on "the unreasonable effectiveness of mathematics" [25] for modeling processes in the real world, and even in disciplines less precise than physics, we often can assume away a surprising amount of detail and still obtain useful answers.

In the remainder of this paper, the mapping of real-world situations to mathematical models is dealt with only in a cursory way. The primary focus is on the purely mathematical aspect of the models, and how Mathematica can be used to solve them.

Relibility modeling

A more complex notion is the hazard function h(t), which is the failure rate, conditional on the unit having survived until a certain point in time. This is also called the "instantaneous age-specific failure rate" or simply the "failure rate".

This conditional rate is based on the general concept of conditional probability: the (unconditional) probability that events A and B both occur is the product of the probability that A occurs, by the probability that B occurs, given that A has occurred (the conditional probability of B given A). In symbols:
    P(A⋂B) = P(A) · P(B|A).
Thus the conditional probability is expressed as

    P(B|A) = .                    (1)

To derive the hazard function, consider the probability that failure occurs between times t and t + Δt:

    P(t <= T <= t + Δt) = f(τ)dτ = F( t + Δt) - F(t) = (1 - F(t)) - (1 - F( t + Δt)) = R( t + Δt) - R(t).
    
Using (1), conditioning the probability on the fact that failure has not occurred before t, i.e., that the unit is still a survivor at time t,

     P(t <= T <= t + Δt | T > t) = =

Note that the hazard function is not a PDF, since ~= 1; it can be shown that R(t)->0 fast enough as t->∞ so that is unbounded

Empirical models of reliability

The blue line is the empirical CDF, which has a step at each ; the step is from (i-1)/n to i/n. The red line is a continuous CDF fitted to the empirical CDF; ways of doing the fitting are described below.

The partition 0 <= 1/n <= . . . i/n <= 1 is biased, since it assumes that P(t < ) = 0, and that is the maximum value for t. There are various methods of compensating for this bias; the preferred method for asymmetrical distributions (where the mean and median are not equal) is Benard's formula:  

      (t <= ) = () = .
    
This is the best approximation for distributions commonly encountered in reliability work. (For a comparison of Benard with other methods in the Weibull case, see [1].)

As an example, take the data set leukemiaRemission, which contains times (in months) to spontaneous remission of leukemia in a control group for a drug study (from [13]). Spontaneous remission of leukemia is apparently common; "failure" here is thus a failure of the carcinoma to continue its uncontrolled growth. The form of the dataset is {. . ., {i, }, . . .}.

Before operating on a set of observation data, taking the data as a discrete function of t, it must be single-valued. The leukemiaRemission dataset is not--some of the remission times occur more than once:

This means that in the kind of step function shown in Figure 1, there is a step greater than one at certain values--e.g.,the step from t = 5 to t = 8 is from 9 to 13. The data is adjusted for this by aggregating the function pairs {10,8}, {11,8}, {12,8}, and {13,8}, leaving only {13,8}.  In general,the aggregation leaves only the observation with the largest i, for a given t. It also swaps the parameters so each pair is {t, i}:

So the dataset,when aggregated, represents a single-valued function from t to the order statistics 1, 2, . . . n:

The aggregated data is used in the following Mathematica function for (t), the empirical CDF, which takes as argument the list of failure time observations. (Note that n is still the number of unaggregated observations.)

The empirical CDF table is always returned. If Plot->True, it also plots the CDF table; PlotWith->{, . . .}, gives a list of other functions to be plotted, for example if a fitted function is to be compared with the empirical CDF (see Figure 3).

Having constructed the empirical CDF, it can be used directly to generate sample failure times; this is illustrated below in the section below on  Simulation or Monte Carlo Methods.  Here we look at deriving other functions from the empirical data--particularly the empirical hazard function,  (t), which is helpful in selecting a distribution type for fitting to the empirical data.

To obtain (t), we start with (t), the estimated PDF:

     () = = = =

The reliability function is () = 1 - (), so from (2),

    () = = = =
    
This last formula can be applied directly to the observation data to derive an empirical (t):

Applying this to leukemiaRemission:

                Figure 2.

The next step is to use this information to determine a distribution for fitting the data.

If the hazard rate is found to be constant, h(t)  = = λ , then R'(t) = - λR(t) , thus R(t) = . Given the constraint R(0) = 1, this solution is unique, so the failure distribution is exponential. In Mathematica terms:

The theoretical basis for the exponential distribution is that the unit has no "memory" of its prior life; the conditional reliability function, given that a unit has survived until time τ, is the same as the reliability function starting at time 0:

    P(T > τ + t | T  > τ) = = =   = R(t) = P(T > t).
    
Empirically, this applies to units such as electronic components that do not wear out. Does it also apply to the leukemia situation? The empirical hazard rate in Figure 2 does not make this clear. Given that it does not seem to be definitely increasing or decreasing, we can look at an average of the first half of the hazard table:

And the last half, with the last observation dropped because it seems like an outlier:

These are not very close, but perhaps (in practical work) close enough to warrant trying a fit to an exponential distribution. First, since () = 1 - (),

Given R(t) = , log(R(t)) = - λt is a straight line. Mathematica's linear regression can be used to find the best fit to the linearized data (IncludeConstant->False is used because we know that there is no constant term):

The value of .96 indicates a good fit to the data, which can be verified visually:

                    Figure 3.

This can now be used as a model for prediction. For example, suppose we are running a trial on a drug that is claimed to induce remission of leukemia. How long should the trial run? Given that spontaneous remissions occur, we may not want the trial to be so long that most patients undergo remission regardless of the drug regimen. Suppose we want a length such that no more than 50% of patients will undergo spontaneous remission. This can be determined by solving R(t) = 0.5. Using the fitted model yields a value of 5.83, so it would be appropriate to run the trial for six months.

Given a continuous distribution with PDF f, or a discrete distribution with probabilities p, the rth moment about the mean μ is or  . The first moment is the mean, the second is the variance, the third is the skewness (a measure of the asymmetry of the distribution), and the fourth is the kurtosis (a measure of how sharply the distribution peaks). The entire set of moments completely characterizes the distribution, and ratios among them (particularly the first four) can be used to test an empirical distribution against a theoretical model.

For example, the coefficient of variation is the standard deviation σ (square root of the variance or second moment) divided by the mean μ; for an exponential distribution this ratio is always 1. Calculating the ratio using the discrete leukemiaRemission data gives

which suggests that an exponential distribution may not be the right model for this data. The coefficient of skewness, defined as the mean of , is always 2 for an exponential distribution, but is less than 1 for the  leukemiaRemission data:  

which again suggests that an exponential distribution may not be the right model for this data, even though the fit is visually good. A better decision could be made by obtaining more data, or by fitting to other types of distribution and picking the best.

Given a set of observations, and an n-parameter distribution hypothesized to be the source of the observations, the maximum likelihood method attempts to find the set of parameter values that maximizes the probability of the given observations. This will be illustrated using the (one-parameter) exponential distribution.

Suppose ,  , . . .,   are believed to be observations from an exponential distribution with parameter λ. Then for each  , the probability density at that value is f() = λ.  Since the observations are assumed to be independent, joint probability is the product of individual probabilities, and the joint density for all the t's is the likelihood function
    L(λ, ,  , . . .,  )  =   .
The value of λ which maximizes L for the given set of t's is the maximum likelihood (ML) estimate.

Since log is monotomically increasing, L has a maximum iff log(L) is maximized. The log likelihood is easier to work with for most distributions used in reliability, so we use
    log(L(λ, ,  , . . .,  ))  =   = n log(λ) - .
    log(L(λ, ,  , . . .,  ))  = -  
                      =  0 when
                          λ  =  n / .
Note that
    log(L(λ, ,  , . . .,  ))  = - <  0
so this is a maximum for the given fixed  ,  , . . .,  .

Applying this to the example above, the maximum likelihood estimate of the exponential rate parameter for the leukemia remission data is:

This can be compared with the linear regression estimate using the Kolmogorov-Smirnov test for goodness of fit (described below):

The smaller K-S statistic for the maximum likelihood estimate indicates it is a better fit. This is in accord with the theory, which says that ML is the least biased estimate of distribution parameters. Nevertheless, in practical reliability work, linear fitting is used more often than ML. ML for the exponential distribution is exceptionally simple; for multiparameter distributions (such as the Weibull, described below) it is more difficult, and apparently practicioners feel that the somewhat better fit does not justify the effort.  In this example, either fit might appear "good enough" visually:   

(The lower continuous curve is the ML estimated fit.)

The Weibull distribution is a flexible two-parameter distribution that is widely used in reliability analysis. The theoretical justification is a theorem that says that in the limit, the smallest value in each of a large number of samples from any distribution will be Weibull-distributed. This weakest link property is valid for many types of units, particularly mechanical parts, where a failure in the weakest part of the unit propagates and causes the entire unit to fail.

The Weibull equations are:

(The Weibull, as well as the exponential distribution, is also available in the Mathematica package Statistics`ContinuousDistributions`.)

In these equations, a is the scale parameter and b is the shape parameter. Note that if b = 1, the result is an exponential distribution with λ = 1/a. Depending on the value of b, a Weibull distribution can have decreasing, constant (exponential), or  increasing hazard rate:

In looking at the hazard rate for a complex system, we can use Weibull distributions and a competing risk model to derive the hazard rate curve for the total system life. The idea is that there is a "burn-in" period with high, but decreasing, hazard rate, as weak components fail and are replaced. Then, during the operating life of the system, the hazard rate is constant; finally, as components mechanically wear out, the hazard rate increases, and presumably the system is retired. Here is a simple competing risk model:

This yields an example of the well-known "bathtub curve" of failure rate over the life of a system.

The ballBearingFailures dataset (from [14]) is failure times (measured in millions of revolutions) of a sample of 23 ball bearings, tested until all had failed. The empirical hazard rate looks like this:   

Smoothing the data makes it easier to interpret:

Aside from one outlying point, the hazard rate seems to increase and then decrease, suggesting more than one failure mode is at work. But nevertheless, we can see how good a fit we get to a single Weibull distribution. Since
    R(t) = exp[-]
    log[R(t)] = -
    log{-log[R(t)]} = b(log t) - b(log a),
letting x = log t, y = log{-log[R(t)]}, we have a linear equation with slope b and intercept -b(log a). After fitting an equation to the data, we can then solve for a and b.

First, we transform the data appropriately, then do the fit:

of .97 indicates a good fit. Note, however, that most empirical datasets don't tell much about the tails of the distribution, very short and very long failure times; and it is the tails that often are needed to differentiate between distributions. E.g., in [26] it is shown that this particular dataset fits either Weibull or Lognormal distributions.

The regression equation yields b = 2.17951, and since intercept = -9.59963 = -b(log a), for a we have:

Visually, this gives a good fit to the data:

As an example of how this model might be used, suppose the manufacturer of the ball bearings in question wishes to determine an economic warranty period. A financial analysis might suggest that no more than 1% of the bearings should need replacement under warranty.  This means we need to find the value of t such that the Weibull reliability function yields R(t) >= 0.99.

This implies that the manufacturer should warrant the bearings for no more than about 10 million revolutions. Given that warranties are ordinarily stated in months, years or miles, some assumption about the duty cycle is needed. Suppose the bearings are used in the hubs of 26-inch radius bicycle wheels. Such a wheel makes about 388 revolutions per mile, so conservatively, the bearings could be warrantied for 25,000 miles. This would be a reasonable life for a bicycle bearing; it would not be very good if the bearings were used in automobile wheels.

Both the leukemia remission and ball bearing datasets were uncensored--that is, all items in the sample were tested to failure. Often, we have censored data of various kinds. In leukemia tests, study participants may drop out before undergoing remission, e.g., because they moved away or died from an unrelated cause. A ball bearing test would typically be run for a fixed length of time or a fixed number of failures, rather than until all items had failed. Thus we would end up with a set of failure times, and another set of times that simply said the items had not failed as of when the test was terminated. Methods exist to deal with censored datasets, but they are beyond the scope of this paper (see [13], [20], [26]).

For the rationale behind computing D+ and D-, see any detailed treatment of the K-S test, e.g., in [12]. The distribution (and thus the critical values) of the K-S statistic depend on the source of the theoretical distribution being tested. If it is fitted based on the empirical data, as is done above, the distribution of the K-S statistic is specific to the type of distibution, e.g., Weibull or exponential.

For the ball bearing data,

Using the K-S table for the Weibull distribution in [12], we find that the null hypothesis can be accepted at a significance level of 0.10; i.e., we don't reject the hypothesis that the data are sampled from the given Weibull distribution. Note that for acceptance of the null hypothesis, a higher significance level is better.

Now let us automate the entire procedure of fitting to a Weibull distribution and testing it for goodness of fit. First, the function fitWeibull fits (and returns) a Weibull CDF from its input, a set or points representing an empirical reliability function. Mathematica's Fit is used here instead of Regress, because we only need the fitted coefficients and not the correlation information. The function returns a list with the Weibull shape and scale parameters, and the CDF as a pure function:

Next, using a table of critical values for the K-S statistic for a fitted Weibull distribution from [12], the goodness of fit is tested directly from the failure dataset. The ksWeibullGof[] function returns the raw data in a list: sample size, Weibull shape and scale parameters, a boolean reject indicating whether the null hypothesis is to be rejected, the significance level alpha, the value of the K-S statistic, and the fitted Weibull CDF as a pure function:

ksWeibullGof[(data_)?ListQ] := 
Module[
{ksWeibullTable, sampleSize, sampleSizeList, fit, cdf, i,shapeParm, scaleParm, empCdf, ksStat, sigList,reject,alpha},
(* Table rows are sample size followed by pairs of 1-alpha, K-S value *) 
ksWeibullTable := {
  {10, {0.9, 0.76}, {0.95, 0.819}, {0.99, 0.944}}, 
  {20, {0.9, 0.779}, {0.95, 0.843}, {0.99, 0.973}}, 
  {50, {0.9, 0.79}, {0.95, 0.856}, {0.99, 0.988}}, 
  {100, {0.9, 0.803}, {0.95, 0.874}, {0.99, 1.007}}}; 
  sampleSize := Length[data];
  sampleSizeList := ksWeibullTable[[Length[ksWeibullTable]]]; 
  empCdf := empiricalCDF[data, Plot -> False]; 
  empCdfSize := Length[empCdf]; 
  fit := fitWeibull[empiricalReliability[data]]; 
  shapeParm := fit[[1]];
  scaleParm := fit[[2]]; 
  cdf := fit[[3]]; 
  (* Select appropriate row from K-S table for sample size *)
  For[
    i = Length[ksWeibullTable], i > 1, i--,
    If[
      sampleSize < ksWeibullTable[[i]][[1]], 
      sampleSizeList = ksWeibullTable[[i - 1]]
    ]
  ]; 
  ksStat := ksStatistic[empCdf, cdf]; 
  adjKsStat := Sqrt[sampleSize]*ksStat; 
  sigList = Null; 
  (* Select nearest higher K-S statistic value from table *)
  For[
    i = 2, i <= Length[sampleSizeList], i++, 
    If[
      adjKsStat < sampleSizeList[[i]][[2]], 
      sigList:= sampleSizeList[[i]]; Break[]
    ]
  ]; 
    reject:=sigList===Null; 
    If[!reject, alpha:=1-sigList[[1]]];
    {sampleSize,shapeParm,scaleParm,reject,alpha,ksStat,cdf}
]

The printKsWeibullGof[] function prints the goodness of fit information,and returns the fitted Weibull function:

An alternative for testing goodness of fit is the Pearson test. The test operates on grouped frequency data, which makes it popular in certain fields, such as social sciences and medicine, where observations are often naturally grouped. It can be cumbersome for data that are not grouped naturally, and is also less satisfactory than Kolmogorov-Smirnov in some contexts for technical reasons (see [2]); in particular, it may not be accurate for distributions with long tails, such as the exponential and Weibull.

The idea is that for grouped data, we compute a statistic that compares , the number of observations found in each group, with , the number we would expect to find, under the hypothesis that the data come from a certain distribution. The statistic is

      =   .

This statistic is approximately -distributed, hence the name of the test. The smaller the value of the statistic, the better the fit. The distribution has one parameter, the number of degrees of freedom. This is the number of independent observations, which for grouped data is the number of groups minus 1.

If the test is of observed data against a hypothesized distribution, is calculated by looking at the largest and smallest observations in the group, say (if we are dealing with failure times) and . If the total number of observations in the data set is n, the expected number in [, ] is n |F() - F()|, where F is the distribution CDF. The null hypothesis in this case is that the observed data come from the hypothesized distribution, and deviations are by chance. The null hypothesis is rejected at the k% level if the probability of the deviation (expressed by the statistic) is less than k%.

Testing the sample failure data involves two steps--grouping the data, and then computing the statistic based on the grouped data and the hypothetical distribution. To achieve acceptable accuracy, the consensus is that there should be at least three groups, with at least five observations per group. Here, any extra observations are added to the last group :

Here the leukemia data is tested against the fitted exponential distribution:

The probability can be evaluated using a table, or using Mathematica's function:

This means that the probability of a value this large or larger is about 32%. Thus if we are (as is typically done) testing at the 5% significance level, we do not reject the null hypothesis. A similar procedure can be carried out on the ball bearing data:

Again, we do not reject the null hypothesis. Note that these are "p values"; the significance level is 1 - p. Note also that since the distribution is known, an exact p-value or significance can be given. In the case of the Kolmogrov-Smirnov test where the hypothetical distribution is fitted from the data being tested, there is no closed form expression or other simple means of computing the distribution, so tables derived by Monte Carlo methods are used. Thus the result is given in terms of critical levels such as 0.10 (or the equivalent p-values such as 0.90) rather than exact values.

There are many systems in the world that exhibit (at least approximately) the following characteristics:

1. At any given point in time, the system is in some defined state, which may be labeled with an integer (if the total number of states is finite or countable) or a real number (if the set of states is continuous).
2. The system periodically undergoes a transition from its current state to another state; in general, the state to which a transition is made is known only in a probabilistic sense. Transitions may occur at finite time intervals, or continuously.
3. The probability of making a transition from state i to state j at time t depends only on i and j; i.e., the system has no memory of its past states. (In the case of a process where transitions may occur continuously in time, we deal with the probability density or rate of transitions, rather than the transition probability.)

These conditions define a Markov process, named after the Russian mathematician who first studied them. Such processes are widely used in applications ranging from modeling traffic queueing, to studying phoneme sequences in natural language. Only a tiny fragment will be dealt with here.

In reliability work, we deal almost exclusively with discrete-state processes, where the states might be whether a unit is operational or not, how many units have failed, the number of repair actions, etc. Typically, the same system can be mapped to more than one state space, and we pick the one that is most useful for the problem at hand. Time is taken to be continuous, since failure can occur at any instant in time. Continuous-time systems are initially defined in terms of discrete time, where transitions occur at intervals Δt;  we can then examine how the process behaves as Δt->0.

To begin with, take a simple example using discrete time. Given a computer system whose availability is only measured to some finite time unit, say an hour, we can build a Markov model in which transitions occur hourly. State 1 is "the computer is up (operational)"; state 2 is "the computer is down (failed)". Suppose that from measurements, the probability that the computer will fail in any given hour is 0.01, and, given that it is down, the probability is .30 that it will be repaired in any given hour after the failure.

We first define the (square) transition probability matrix P, such that is the probability of making a transition from state i to state j. For this case,

Note the use of the constraint that the sum of entries in each column vector must be 1, since the set of transitions from state i is mutually exclusive and exhaustive. Now, we add a state vector x giving the probabilities of being in the various states at a point in time (again, the entries must sum to 1); initially, we can assume that the system is up, so x(0) = (1, 0). Since is the probability of being in state i before the transition, the probability of being in state j after the transition is  . This is the jth entry in the vector Px. For example, after the first transition in the example, the state vector is

After two transitions, it is

In general, after n transitions the state vector is P(P(P...x)...) = x. In situations such as this example, an obvious question is whether the process is stationary, i.e., whether the state vector converges to some value over time. Given a tool like Mathematica, the "quick and dirty" way to find out is to determine the stationary value empirically, by computing powers of P. If x is taken as a vector in which the ith entry is 1 and all others are 0,  x will be the ith column of  . In this case, P converges:

Thus over a long period of time, we expect the system to be available about 97% of the time. In "normal" systems encountered in reliability the process is always stationary, and it can be shown that, regardless of the initial state, the state vector converges to the (identical) column vectors of  , as n->∞. (See, e.g., [10], [19]).

A more formal method of determining the value of the column vectors is shown below, but first it is interesting to look at a closed-form solution for this example, using differential equations.

From the conditions of a Markov process, the probability (or probability density) of transitions does not change with time. Making the example more realistic, instead of saying that the probability of transition is, say, λ on the hour, we say the rate of transition is λ per hour. Since the rate is constant, for any interval Δt the probability of transition within the interval is λ Δt. Then the transition matrix (using λ for the rate of failure and μ for the rate of repair), for transitions every Δt, becomes   

    P = ,  thus for one transition the equation for the state vector is
    
    x(t+Δt)  =   =  Px(t)  = .

Separating the two components of the state vector, and after some algebra, this reduces to

     = - λ (t) + μ (t)     and     = λ (t) - μ (t).

Taking limits as Δt -> 0,

     ' (t) = - λ (t) + μ (t),    ' (t) = λ (t) - μ (t).

Finally, taking advantage of the fact that (t) = 1 - (t), after some algebra,

     ' (t) =  - (λ + μ) (t) + μ ,         ' (t) =  - (λ + μ) (t) + λ .
     
These can be solved by Mathematica, using the example values λ = 0.01, μ = 0.30, and the initial conditions (0) = 1, (0) = 0:  

This is identical to the solution computed above using matrix powers.

Now consider a more complex process, shown in figure 4 below. There are two redundant systems, A and B; either can fail independently, and if down, can be repaired and brought back up. The figure shows the relevant transition rates (the probability of a direct transition from state 1 to state 4 is small enough to be neglected). The most important question is the probability of total system outage, represented by state 4.

                                Figure 4.

The transition matrix is

As in the previous example, we can find the stationary probabilities by brute force:

(In this case, Mathematica does not seem to be able to handle  Limit[MatrixPower[case2P,n],n->∞]--the result is a zero matrix. It does converge to the value above, however, as can be verified by computing a series of large powers of the matrix.)

More formally, in finding the stationary probabilities, we are looking for the state vector x such that Px = x. In other words, an eigenvector x corresponding to the eigenvalue 1, and satisfying   = 1, since the are probabilities of a set of mutually exclusive and exhaustive states. Starting with

The desired eigenvector is the first:

However, the components do not sum to one. Since an eigenvector is unuque only up to a scalar multiple, it can be normalized by the sum of its components to produce the desired result:

We can confirm that it is the desired eigenvector:

And also note that it is identical to the result produced above by matrix powers.

Markov models are limited by the fact that the "memoryless" property implies exponential transition rates; thus modeling of, say, Weibull failure rates is ruled out. But during the operating life of many systems (the flat part of the bathtub curve), exponential rates are a good enough approximation to make the simplicity of the Markov approach appealing.

Simulation or Monte Carlo methods

One aspect of simulation modeling is defining the network model, which is specific to the situation under study. A more general problem, arising in all simulation models, is this: given a random variable X characterized by a probability distribution, to generate a time series of values , , , . . . from this distribution. All languages used for writing simulations, including Mathematica, supply methods for generating random numbers from a uniform distribution. By definition, if U is a uniform random variable on the interval [a, b], the probability of a value of U being found in the interval [u, u + Δu] is . Recalling the discussion above related to figure 1, for any random variable X, the probability P that the random variable has a value <= x is uniformly distributed on the interval [0, 1]. So if the CDF of X has an inverse : P -> X, we can use the inverse function to map the uniform variable U[0, 1] to X[0, ∞].

For example, take the CDF for the exponential distribution,  p = F(t) = 1 - .  The inverse function is t = , but note that p and 1-p have the same distribution.  So this function will return a random exponential variate for the parameter λ :

Here the function is used to generate a data set simulating the leukemia remission data, which is plotted against the continuous CDF:

(Compare this to the original data plot.)

This might be used to set up a simulation model that showed how natural remissions competed with the hypothesized effect of a new drug, in order to validate the design of a trial.

Since a simulation might need to generate thousands of values from a given distribution, computing the inverse function for each value may impact performance, particularly for CDFs (such as the normal) which do not have closed-form inverses. Given that the results only have a finite precision, standard practice is to initialize a table mapping p to t in increments, say of .01. This is similar to the empirical CDF, and in fact an empirical CDF could be used in a simulation model for this purpose.

Taking as an example the Weibull distribution for the ball bearing data, here is a table-driven generator of Weibull random variates. First a table is built. In order to insure coverage of the significant portion of the probability space, the function first computes the value of t with a cumulative probability of 0.999; then it iterates from 0 through this t, computing at least 100 values (the number of values could be increased for more accuracy). Then the table is used to generate the random variates, by mapping a value from U[0,1] to the closest table value.

Here a table is created to simulate the ball bearing failure data

and a simulated data set is created from the table and plotted against the CDF:

(Compare this to the original data plot.)

To illustrate a simple, but more complete, simulation model, take the first example from the section on Markov models. As shown in the figure below, a system may fail, and having failed, it may be repaired. The failure time and repair time distributions are exponential, with the indicated rate parameters:

The general scheme of the model is, starting in state 1, pick a random failure time. All time prior to the failure was spent in state 1, so add this time to the cumulative time in state 1, and move to state 2. Then pick a random repair time, add it to the cumulative time in state 2, and move to state 1. This process is iterated some number of times, and then results are computed; for example,    gives the long-term probability that the system is operational. Here is the model function, which takes as parameters the failure and repair rates, number of iterations, and initial state:

And here it is run with the rate parameters from the example:

Comparing this to the analytic solution {0.967742, 0.0322581}, derived above,

shows it is accurate to within 1%. Convergence to the expected value is quite slow as the number of iterations is increased.

Methods exist for computing statistical confidence intervals on the outputs, but these become less tractable as the model becomes realistically large. Typically, engineers are content to perform a sensitivity analysis by varying parameters over a range they consider realistic, and observing how the variations affect the output.

Summary

References

[1] Abernathy, Bob, and Wes Fulton. Weibull News 11, 1997. http://www.weibullnews.com/wnews11.html.

[3] Department of Defense. MIL-HDBK-217F: Reliability Prediction of Electronic Equipment. Washington, D. C.: Department of Defense, 1991. (Downloadable from http://astimage.daps.dla.mil/quicksearch/.)

[4] Department of Defense. MIL-HDBK-338B: Electronic Reliability Design Handbook. Washington, D. C.: Department of Defense, 1998. (Downloadable from http://astimage.daps.dla.mil/quicksearch/.)

[5]  Feller, William. An Introduction to Probability Theory and Its Applications, Volume I. New York, N. Y.: John Wiley & Sons, 1968.

[8]  InfoPlease. "United States - U.S. Statistics - Population". Available online at http://www.infoplease.com/ipa/A0004997.html.

[9]  Inter-university Consortium for Political and Social Research. "Historical Demographic, Economic and Social Data: The United States, 1790-1970". Ann Arbor, Mich.: n.d. Data extract available online at http://fisher.lib.virginia.edu/census/.

[15]  McLinn, James A. Weibull Analysis Primer. American Society for Quality Control Reliability Division monograph, 1997.

[16]  Meridian Marketing Group. "Get top performance from your system with simulation" (Advertisement for GPSS/H).
http://www.meridian-marketing.com/GPSS_H/wo_p1.html.

[17]  Musa, John D., et al. Software Reliability: Measurement, Prediction, Application. New York, N. Y.: McGraw-Hill, 1987.

[21]  Pukite, Jan, and Paul Pukite. Modeling for Reliability Analysis: Markov Modeling for Reliability, Maintainability, Safety, and Supportability Analysis of Complex Computer Systems. New York, N. Y.: IEEE Press, 1998.

[22]  Sobol', Ilya M. A Primer for the Monte Carlo Method. Boca Raton, Fla.: CRC Press, 1994.

[24]  U. S. Census Bureau. "First Census 2000 Results". Available online at http://www.census.gov/main/www/cen2000.html.