Parameter Estimation

Parameter estimation refers to the process of using sample data (in our case times-to-failure or success data) to estimate the parameters of the selected distribution. Several parameter estimation methods are available. This section presents an overview of the available methods used in life data analysis. More specifically, we start with the relatively simple method of probability plotting and continue with the more sophisticated methods of Rank Regression (or Least Squares) and Maximum Likelihood.

Probability Plotting

The least mathematically intensive method for parameter estimation is the method of probability plotting. As the term implies, probability plotting involves a physical plot of the data on specially constructed probability plotting paper. This method is easily implemented by hand, given that one can obtain the appropriate probability plotting paper.

Illustrating the Method for the 2-Parameter Weibull Distribution

The method of probability plotting takes the cdf of the distribution and attempts to linearize it by employing a specially constructed paper. This is best illustrated using the 2-parameter Weibull distribution.

In the case of the two-parameter Weibull distribution, the cdf (also the unreliability Q(t)) is given by:

[math]\displaystyle{ F(t)=Q(t)=1-{e^{-\left(\tfrac{t}{\eta}\right)^{\beta}}} }[/math]

Linearizing the Weibull Unreliability Function
This function can then be linearized (i.e., put in the common form of y = m'x + b format) as follows:

[math]\displaystyle{ \begin{align} Q(t)= & 1-{e^{-\left(\tfrac{t}{\eta}\right)^{\beta}}} \\ \ln (1-Q(t))= & \ln \left[ {e^{-\left(\tfrac{t}{\eta}\right)^{\beta}}} \right] \\ \ln (1-Q(t))=& -\left(\tfrac{t}{\eta}\right)^{\beta} \\ \ln ( -\ln (1-Q(t)))= & \beta \left(\ln \left( \frac{t}{\eta }\right)\right) \\ \ln \left( \ln \left( \frac{1}{1-Q(t)}\right) \right) = & \beta\ln{ t} -\beta(\eta ) \\ \end{align} }[/math]

Then by setting

[math]\displaystyle{ y=\ln \left( \ln \left( \frac{1}{1-Q(t)} \right) \right) }[/math]

and

[math]\displaystyle{ x=\ln \left( t \right) }[/math]

the equation can then be rewritten as:

[math]\displaystyle{ y=\beta x-\beta \ln \left( \eta \right) }[/math]

which is now a linear equation with a slope of

Slope = m = β

and an intercept of

Intercept[math]\displaystyle{ =b=-\beta \cdot ln(\eta) }[/math].

Constructing the Paper

The next task is to construct the Weibull probability plotting paper with the appropriate y - and x-axes. The x-axis trasnformation is simply logarithmic. The y-axis, is a bit more complex requiring a double log reciprocal transformation, or 

[math]\displaystyle{ y=\ln \left(\ln \left( \frac{1}{1-Q(t)} ) \right) \right) }[/math]

where Q(t) is the unreliability.

Such papers have been created by different vendors and are called probability plotting papers. Weibull.com has different plotting papers available for download.

To illustrate, consider the following probability plot on a slightly different type of Weibull probability paper.

This paper is constructed based on the mentioned y - and x - transformations, where the y-axis represents unreliability and the x-axis represents time. Both of these values must be known for each time-to-failure point we want to plot.

Then, given the y and x value for each point, the points can easily be put on the plot. Once the points have been placed on the plot, the best possible straight line is drawn through these points. Once the line has been drawn, the slope of the line can be obtained (some probability papers include a slope indicator to simplify this calculation). This is the parameter β, which is the value of the slope. To determine the scale parameter, η (also called the characteristic life), one reads the time from the x-axis corresponding to Q(t)=63.2%.

Note that at

[math]\displaystyle{ \begin{align} Q(t=\eta)= & 1-{{e}^{-{{\left( \tfrac{t}{\eta } \right)}^{\beta }}}} \\ = & 1-{{e}^{-1}} \\ = & 0.632 \\ = & 63.2% \end{align} }[/math]

Thus, if we enter the y axis at Q(t)=63.2%, the corresponding value of t will be equal to η. Thus, using this simple methodology, the parameters of the Weibull distribution can be estimated.

Determining the x and y Position of the Plot Points

The points on the plot represent our data or, more specifically, our times-to-failure data. If, for example, we tested four units that failed at 10, 20, 30 and 40 hours, then we would use these times as our x values or time values.

Determining the appropriate y plotting positions, or the unreliability values, is a little more complex. To determine the y plotting positions, we must first determine a value indicating the corresponding unreliability for that failure. In other words, we need to obtain the cumulative percent failed for each time-to-failure. For example, the cumulative percent failed by 10 hours may be 25%, by 20 hours 50%, and so forth. This is a simple method illustrating the idea. The problem with this simple method is the fact that the 100% point is not defined on most probability plots; thus, an alternative and more robust approach must be used. The most widely used method of determining this value is the method of obtaining the median rank for each failure. This is discussed next.

Median Ranks

Median ranks are used to obtain an estimate of the unreliability for each failure. It is the value that the true probability of failure, Q(T_j), should have at the j^t'h failure out of a sample of N units at the 50% confidence level.

The rank can be found for any percentage point, P, greater than zero and less than one, by solving the cumulative binomial equation for Z . This represents the rank, or unreliability estimate, for the j^t'h failure in the following equation for the cumulative binomial:

[math]\displaystyle{ P=\underset{k=j}{\overset{N}{\mathop \sum }}\,\left( \begin{matrix} N \\ k \\ \end{matrix} \right){{Z}^{k}}{{\left( 1-Z \right)}^{N-k}} }[/math]

where N is the sample size and j the order number.

The median rank is obtained by solving this equation for Z at P = 0.50,

[math]\displaystyle{ 0.50=\underset{k=j}{\overset{N}{\mathop \sum }}\,\left( \begin{matrix} N \\ k \\ \end{matrix} \right){{Z}^{k}}{{\left( 1-Z \right)}^{N-k}} }[/math]

For example, if N=4 and we have four failures, we would solve the median rank equation for the value of Z four times; once for each failure with j= 1, 2, 3 and 4. This result can then be used as the unreliability estimate for each failure or the y plotting position. (See also the chapter The Weibull distribution for a step-by-step example of this method.) The solution of cumuative binomial equation for Z requires the use of numerical methods.

Beta and F Distributions Approach

A more straightforward and easier method of estimating median ranks is by applying two transformations to the cumulative binomial equation, first to the beta distribution and then to the F distribution, resulting in [12, 13],

[math]\displaystyle{ \begin{array}{*{35}{l}} MR & = & \tfrac{1}{1+\tfrac{N-j+1}{j}{{F}_{0.50;m;n}}} \\ m & = & 2(N-j+1) \\ n & = & 2j \\ \end{array} }[/math]

where F_0.50;m;n denotes the F distribution at the 0.50 point, with m and n degrees of freedom, for failure j out of N units.

Benard's Approximation for Median Ranks

Another quick, and less accurate, approximation of the median ranks is also given by:

[math]\displaystyle{ MR = \frac{{j - 0.3}}{{N + 0.4}} }[/math]

This approximation of the median ranks is also known as Benard's approximation.

Kaplan-Meier

The Kaplan-Meier estimator (also known as the product limit estimator) is used as an alternative to the median ranks method for calculating the estimates of the unreliability for probability plotting purposes. The equation of the estimator is given by:

[math]\displaystyle{ \widehat{F}({{t}_{i}})=1-\underset{j=1}{\overset{i}{\mathop \prod }}\,\frac{{{n}_{j}}-{{r}_{j}}}{{{n}_{j}}},\text{ }i=1,...,m }[/math]

where:

[math]\displaystyle{ \begin{align} m = & {\text{total number of data points}} \\ n = & {\text{the total number of units}} \\ {n_i} = & n - \sum_{j = 0}^{i - 1}{s_j} - \sum_{j = 0}^{i - 1}{r_j}, \text{i = 1,...,m }\\ {r_j} = & {\text{ number of failures in the }}{j^{th}}{\text{ data group, and}} \\ {s_j} = & {\text{number of surviving units in the }}{j^{th}}{\text{ data group}} \\ \end{align} }[/math]

Analysis of Right Censored (Suspended) Data

All available data should be considered in the analysis of times-to-failure data. This includes the case when a particular unit in a sample has been removed from the test prior to failure. An item, or unit, which is removed from a reliability test prior to failure, or a unit which is in the field and is still operating at the time the reliability of these units is to be determined, is called a suspended item or right censored observation or right censored data point. Suspended items analysis would also be considered when:

1. We need to make an analysis of the available results before test completion.
2. The failure modes which are occurring are different than those anticipated and such units are withdrawn from the test.
3. We need to analyze a single mode and the actual data set comprises multiple modes.
4. A warranty analysis is to be made of all units in the field (non-failed and failed units). The non-failed units are considered to be suspended items (or right censored).

Probability Plotting and Rank Regression Analysis of Right Censored Data

When using the probability plotting or rank regression method to accommodate the fact that units in the data set did not fail, or were suspended, we need to adjust their probability of failure, or unreliability. As discussed before, estimates of the unreliability for complete data are obtained using the median ranks approach. The following methodology illustrates how adjusted median ranks are computed to account for right censored data. To better illustrate the methodology, consider the following example in  [20] where five items are tested resulting in three failures and two suspensions.

Item Number(Position)	State*"F" or "S"	Life of item, hr
1	F1	5,100
2	S1	9,500
3	F2	15,000
4	S2	22,000
5	F3	40,000

The methodology for plotting suspended items involves adjusting the rank positions and plotting the data based on new positions, determined by the location of the suspensions. If we consider these five units, the following methodology would be used: The first item must be the first failure; hence, it is assigned failure order number j = 1. The actual failure order number (or position) of the second failure, F₂, is in doubt. It could either be in position 2 or in position 3. Had S₁ not been withdrawn from the test at 9,500 hours, it could have operated successfully past 15,000 hours, thus placing F₂ in position 2. Alternatively, S₁ could also have failed before 15,000 hours, thus placing F₂ in position 3. In this case, the failure order number for F₂ will be some number between 2 and 3. To determine this number, consider the following:

We can find the number of ways the second failure can occur in either order number 2 (position 2) or order number 3 (position 3). The possible ways are listed next.

F2 in Position 2						OR	F2 in Position 3
1	2	3	4	5	6		1	2
F1	F1	F1	F1	F1	F1		F1	F1
F2	F2	F2	F2	F2	F2		S1	S1
S1	S2	F3	S1	S2	F3		F2	F2
S2	S1	S1	F3	F3	S2		S2	F3
F3	F3	S2	S2	S1	S1		F3	S2

It can be seen that F₂ can occur in the second position six ways and in the third position two ways. The most probable position is the average of these possible ways, or the mean order number ( MON ), given by:

[math]\displaystyle{ {{F}_{2}}=MO{{N}_{2}}=\frac{(6\times 2)+(2\times 3)}{6+2}=2.25 }[/math]

Using the same logic on the third failure, it can be located in position numbers 3, 4 and 5 in the possible ways listed next.

F3 in Position 3		OR	F3 in Position 4			OR	F3 in Position 5
1	2		1	2	3		1	2	3
F1	F1		F1	F1	F1		F1	F1	F1
F2	F2		S1	F2	F2		S1	F2	F2
F3	F3		F2	S1	S2		F2	S1	S2
S1	S2		F3	F3	F3		S2	S2	S1
S2	S1		S2	S2	S1		F3	F3	F3

Then, the mean order number for the third failure, (Item 5) is:

[math]\displaystyle{ MO{{N}_{3}}=\frac{(2\times 3)+(3\times 4)+(3\times 5)}{2+3+3}=4.125 }[/math]

Once the mean order number for each failure has been established, we obtain the median rank positions for these failures at their mean order number. Specifically, we obtain the median rank of the order numbers 1, 2.25 and 4.125 out of a sample size of 5, as given next.

Plotting Positions for the Failures(Sample Size=5)
Failure Number	MON	Median Rank Position(%)
1:F1	1	13%
2:F2	2.25	36%
3:F3	4.125	71%

Once the median rank values have been obtained, the probability plotting analysis is identical to that presented before. As you might have noticed, this methodology is rather laborious. Other techniques and shortcuts have been developed over the years to streamline this procedure. For more details on this method, see Kececioglu [20].

Shortfalls of this Rank Adjustment Method

Even though the rank adjustment method is the most widely used method for performing suspended items analysis, we would like to point out the following shortcoming. As you may have noticed from this analysis of suspended items, only the position where the failure occurred is taken into account, and not the exact time-to-suspension. For example, this methodology would yield the exact same results for the next two cases.

Case 1			Case 2
Item Number	State*"F" or "S"	Life of an item, hr	Item number	State*,"F" or "S"	Life of item, hr
1	F1	1,000	1	F1	1,000
2	S1	1,100	2	S1	9,700
3	S2	1,200	3	S2	9,800
4	S3	1,300	4	S3	9,900
5	F2	10,000	5	F2	10,000
* F − Failed, S − Suspended			* F − Failed, S − Suspended

This shortfall is significant when the number of failures is small and the number of suspensions is large and not spread uniformly between failures, as with these data. In cases like this, it is highly recommended that one use maximum likelihood estimation (MLE) to estimate the parameters instead of using least squares, since maximum likelihood does not look at ranks or plotting positions, but rather considers each unique time-to-failure or suspension. For the data given above, the results are as follows. The estimated parameters using the method just described are the same for both cases (1 and 2):

[math]\displaystyle{ \begin{array}{*{35}{l}} \widehat{\beta }= & \text{0}\text{.81} \\ \widehat{\eta }= & \text{11,417 hr} \\ \end{array} }[/math]

However, the MLE results for Case 1 are:

[math]\displaystyle{ \begin{array}{*{35}{l}} \widehat{\beta }= & \text{1}\text{.33} \\ \widehat{\eta }= & \text{6,900 hr} \\ \end{array} }[/math]

and the MLE results for Case 2 are:

[math]\displaystyle{ \begin{array}{*{35}{l}} \widehat{\beta }= & \text{0}\text{.9337} \\ \widehat{\eta }= & \text{21,348 hr} \\ \end{array} }[/math]

As we can see, there is a sizable difference in the results of the two sets calculated using MLE and the results using regression. The results for both cases are identical when using the regression estimation technique, as regression considers only the positions of the suspensions. The MLE results are quite different for the two cases, with the second case having a much larger value of η, which is due to the higher values of the suspension times in Case 2. This is because the maximum likelihood technique, unlike rank regression, considers the values of the suspensions when estimating the parameters. This is illustrated in the section MLE Analysis of Right Censored Data.

Probability Plots for Other Distributions

This same methodology can be applied to other distributions which have cdf equations that can be linearized. Different probability papers exist for each distribution, since different distributions have different cdf equations. Weibull++ automatically creates these plots for you when choosing a probability plot for a particular distribution. Special scales on these plots allow the parameter estimates to be derived directly from the plots, similar to the way β and η were obtained from the Weibull probability plot. These are discussed in the individual distribution chapters.

Let's assume six identical units are reliability tested at the same application and operation stress levels. All of these units fail during the test after operating for the following times (in hours): 96, 257, 498, 763, 1051 and 1744.

The steps for using the probability plotting method to determine the parameters of the exponential pdf representing the data are as follows:

Rank the times-to-failure in ascending order as shown next.

[math]\displaystyle{ \begin{matrix} \text{Failure} & \text{Failure Order Number} \\ \text{Time (Hr)} & \text{out of a Sample Size of 6} \\ \text{96} & \text{1} \\ \text{257} & \text{2} \\ \text{498} & \text{3} \\ \text{763} & \text{4} \\ \text{1,051} & \text{5} \\ \text{1,744} & \text{6} \\ \end{matrix}\,\! }[/math]

Obtain their median rank plotting positions. Median rank positions are used instead of other ranking methods because median ranks are at a specific confidence level (50%).

The times-to-failure, with their corresponding median ranks, are shown next:

[math]\displaystyle{ \begin{matrix} \text{Failure} & \text{Median} \\ \text{Time (Hr)} & \text{Rank, }% \\ \text{96} & \text{10}\text{.91} \\ \text{257} & \text{26}\text{.44} \\ \text{498} & \text{42}\text{.14} \\ \text{763} & \text{57}\text{.86} \\ \text{1,051} & \text{73}\text{.56} \\ \text{1,744} & \text{89}\text{.10} \\ \end{matrix}\,\! }[/math]

On an exponential probability paper, plot the times on the x-axis and their corresponding rank value on the y-axis. The next figure displays an example of an exponential probability paper. The paper is simply a log-linear paper.

Draw the best possible straight line that goes through the [math]\displaystyle{ t=0\,\! }[/math] and [math]\displaystyle{ (t)=100%\,\! }[/math] point and through the plotted points (as shown in the plot below).

At the [math]\displaystyle{ Q(t)=63.2%\,\! }[/math] or [math]\displaystyle{ R(t)=36.8%\,\! }[/math] ordinate point, draw a straight horizontal line until this line intersects the fitted straight line. Draw a vertical line through this intersection until it crosses the abscissa. The value at the intersection of the abscissa is the estimate of the mean. For this case, [math]\displaystyle{ \widehat{\mu }=833\,\! }[/math] hours which means that [math]\displaystyle{ \lambda =\tfrac{1}{\mu }=0.0012\,\! }[/math] (This is always at 63.2% because [math]\displaystyle{ (T)=1-{{e}^{-\tfrac{\mu }{\mu }}}=1-{{e}^{-1}}=0.632=63.2%)\,\! }[/math].

Now any reliability value for any mission time [math]\displaystyle{ t\,\! }[/math] can be obtained. For example, the reliability for a mission of 15 hours, or any other time, can now be obtained either from the plot or analytically.

To obtain the value from the plot, draw a vertical line from the abscissa, at [math]\displaystyle{ t=15\,\! }[/math] hours, to the fitted line. Draw a horizontal line from this intersection to the ordinate and read [math]\displaystyle{ R(t)\,\! }[/math]. In this case, [math]\displaystyle{ R(t=15)=98.15%\,\! }[/math]. This can also be obtained analytically, from the exponential reliability function.

Some Shortfalls of Manual Probability Plotting

Besides the most obvious drawback to probability plotting, which is the amount of effort required, manual probability plotting is not always consistent in the results. Two people plotting a straight line through a set of points will not always draw this line the same way, and thus will come up with slightly different results. This method was used primarily before the widespread use of computers that could easily perform the calculations for more complicated parameter estimation methods, such as the least squares and maximum likelihood methods.

Parameter estimation refers to the process of using sample data (in our case times-to-failure or success data) to estimate the parameters of the selected distribution. Several parameter estimation methods are available. This section presents an overview of the available methods used in life data analysis. More specifically, we start with the relatively simple method of probability plotting and continue with the more sophisticated methods of Rank Regression (or Least Squares) and Maximum Likelihood.

Probability Plotting

The least mathematically intensive method for parameter estimation is the method of probability plotting. As the term implies, probability plotting involves a physical plot of the data on specially constructed probability plotting paper. This method is easily implemented by hand, given that one can obtain the appropriate probability plotting paper.

Illustrating the Method for the 2-Parameter Weibull Distribution

The method of probability plotting takes the cdf of the distribution and attempts to linearize it by employing a specially constructed paper. This is best illustrated using the 2-parameter Weibull distribution.

In the case of the two-parameter Weibull distribution, the cdf (also the unreliability Q(t)) is given by:

[math]\displaystyle{ F(t)=Q(t)=1-{e^{-\left(\tfrac{t}{\eta}\right)^{\beta}}} }[/math]

Linearizing the Weibull Unreliability Function
This function can then be linearized (i.e., put in the common form of y = m'x + b format) as follows:

[math]\displaystyle{ \begin{align} Q(t)= & 1-{e^{-\left(\tfrac{t}{\eta}\right)^{\beta}}} \\ \ln (1-Q(t))= & \ln \left[ {e^{-\left(\tfrac{t}{\eta}\right)^{\beta}}} \right] \\ \ln (1-Q(t))=& -\left(\tfrac{t}{\eta}\right)^{\beta} \\ \ln ( -\ln (1-Q(t)))= & \beta \left(\ln \left( \frac{t}{\eta }\right)\right) \\ \ln \left( \ln \left( \frac{1}{1-Q(t)}\right) \right) = & \beta\ln{ t} -\beta(\eta ) \\ \end{align} }[/math]

Then by setting

[math]\displaystyle{ y=\ln \left( \ln \left( \frac{1}{1-Q(t)} \right) \right) }[/math]

and

[math]\displaystyle{ x=\ln \left( t \right) }[/math]

the equation can then be rewritten as:

[math]\displaystyle{ y=\beta x-\beta \ln \left( \eta \right) }[/math]

which is now a linear equation with a slope of

Slope = m = β

and an intercept of

Intercept[math]\displaystyle{ =b=-\beta \cdot ln(\eta) }[/math].

Constructing the Paper

The next task is to construct the Weibull probability plotting paper with the appropriate y - and x-axes. The x-axis trasnformation is simply logarithmic. The y-axis, is a bit more complex requiring a double log reciprocal transformation, or 

[math]\displaystyle{ y=\ln \left(\ln \left( \frac{1}{1-Q(t)} ) \right) \right) }[/math]

where Q(t) is the unreliability.

Such papers have been created by different vendors and are called probability plotting papers. Weibull.com has different plotting papers available for download.

To illustrate, consider the following probability plot on a slightly different type of Weibull probability paper.

This paper is constructed based on the mentioned y - and x - transformations, where the y-axis represents unreliability and the x-axis represents time. Both of these values must be known for each time-to-failure point we want to plot.

Then, given the y and x value for each point, the points can easily be put on the plot. Once the points have been placed on the plot, the best possible straight line is drawn through these points. Once the line has been drawn, the slope of the line can be obtained (some probability papers include a slope indicator to simplify this calculation). This is the parameter β, which is the value of the slope. To determine the scale parameter, η (also called the characteristic life), one reads the time from the x-axis corresponding to Q(t)=63.2%.

Note that at

[math]\displaystyle{ \begin{align} Q(t=\eta)= & 1-{{e}^{-{{\left( \tfrac{t}{\eta } \right)}^{\beta }}}} \\ = & 1-{{e}^{-1}} \\ = & 0.632 \\ = & 63.2% \end{align} }[/math]

Thus, if we enter the y axis at Q(t)=63.2%, the corresponding value of t will be equal to η. Thus, using this simple methodology, the parameters of the Weibull distribution can be estimated.

Determining the x and y Position of the Plot Points

The points on the plot represent our data or, more specifically, our times-to-failure data. If, for example, we tested four units that failed at 10, 20, 30 and 40 hours, then we would use these times as our x values or time values.

Determining the appropriate y plotting positions, or the unreliability values, is a little more complex. To determine the y plotting positions, we must first determine a value indicating the corresponding unreliability for that failure. In other words, we need to obtain the cumulative percent failed for each time-to-failure. For example, the cumulative percent failed by 10 hours may be 25%, by 20 hours 50%, and so forth. This is a simple method illustrating the idea. The problem with this simple method is the fact that the 100% point is not defined on most probability plots; thus, an alternative and more robust approach must be used. The most widely used method of determining this value is the method of obtaining the median rank for each failure. This is discussed next.

Median Ranks

Median ranks are used to obtain an estimate of the unreliability for each failure. It is the value that the true probability of failure, Q(T_j), should have at the j^t'h failure out of a sample of N units at the 50% confidence level.

The rank can be found for any percentage point, P, greater than zero and less than one, by solving the cumulative binomial equation for Z . This represents the rank, or unreliability estimate, for the j^t'h failure in the following equation for the cumulative binomial:

[math]\displaystyle{ P=\underset{k=j}{\overset{N}{\mathop \sum }}\,\left( \begin{matrix} N \\ k \\ \end{matrix} \right){{Z}^{k}}{{\left( 1-Z \right)}^{N-k}} }[/math]

where N is the sample size and j the order number.

The median rank is obtained by solving this equation for Z at P = 0.50,

[math]\displaystyle{ 0.50=\underset{k=j}{\overset{N}{\mathop \sum }}\,\left( \begin{matrix} N \\ k \\ \end{matrix} \right){{Z}^{k}}{{\left( 1-Z \right)}^{N-k}} }[/math]

For example, if N=4 and we have four failures, we would solve the median rank equation for the value of Z four times; once for each failure with j= 1, 2, 3 and 4. This result can then be used as the unreliability estimate for each failure or the y plotting position. (See also the chapter The Weibull distribution for a step-by-step example of this method.) The solution of cumuative binomial equation for Z requires the use of numerical methods.

Beta and F Distributions Approach

A more straightforward and easier method of estimating median ranks is by applying two transformations to the cumulative binomial equation, first to the beta distribution and then to the F distribution, resulting in [12, 13],

[math]\displaystyle{ \begin{array}{*{35}{l}} MR & = & \tfrac{1}{1+\tfrac{N-j+1}{j}{{F}_{0.50;m;n}}} \\ m & = & 2(N-j+1) \\ n & = & 2j \\ \end{array} }[/math]

where F_0.50;m;n denotes the F distribution at the 0.50 point, with m and n degrees of freedom, for failure j out of N units.

Benard's Approximation for Median Ranks

Another quick, and less accurate, approximation of the median ranks is also given by:

[math]\displaystyle{ MR = \frac{{j - 0.3}}{{N + 0.4}} }[/math]

This approximation of the median ranks is also known as Benard's approximation.

Kaplan-Meier

The Kaplan-Meier estimator (also known as the product limit estimator) is used as an alternative to the median ranks method for calculating the estimates of the unreliability for probability plotting purposes. The equation of the estimator is given by:

[math]\displaystyle{ \widehat{F}({{t}_{i}})=1-\underset{j=1}{\overset{i}{\mathop \prod }}\,\frac{{{n}_{j}}-{{r}_{j}}}{{{n}_{j}}},\text{ }i=1,...,m }[/math]

where:

[math]\displaystyle{ \begin{align} m = & {\text{total number of data points}} \\ n = & {\text{the total number of units}} \\ {n_i} = & n - \sum_{j = 0}^{i - 1}{s_j} - \sum_{j = 0}^{i - 1}{r_j}, \text{i = 1,...,m }\\ {r_j} = & {\text{ number of failures in the }}{j^{th}}{\text{ data group, and}} \\ {s_j} = & {\text{number of surviving units in the }}{j^{th}}{\text{ data group}} \\ \end{align} }[/math]

Analysis of Right Censored (Suspended) Data

All available data should be considered in the analysis of times-to-failure data. This includes the case when a particular unit in a sample has been removed from the test prior to failure. An item, or unit, which is removed from a reliability test prior to failure, or a unit which is in the field and is still operating at the time the reliability of these units is to be determined, is called a suspended item or right censored observation or right censored data point. Suspended items analysis would also be considered when:

1. We need to make an analysis of the available results before test completion.
2. The failure modes which are occurring are different than those anticipated and such units are withdrawn from the test.
3. We need to analyze a single mode and the actual data set comprises multiple modes.
4. A warranty analysis is to be made of all units in the field (non-failed and failed units). The non-failed units are considered to be suspended items (or right censored).

Probability Plotting and Rank Regression Analysis of Right Censored Data

When using the probability plotting or rank regression method to accommodate the fact that units in the data set did not fail, or were suspended, we need to adjust their probability of failure, or unreliability. As discussed before, estimates of the unreliability for complete data are obtained using the median ranks approach. The following methodology illustrates how adjusted median ranks are computed to account for right censored data. To better illustrate the methodology, consider the following example in  [20] where five items are tested resulting in three failures and two suspensions.

Item Number(Position)	State*"F" or "S"	Life of item, hr
1	F1	5,100
2	S1	9,500
3	F2	15,000
4	S2	22,000
5	F3	40,000

The methodology for plotting suspended items involves adjusting the rank positions and plotting the data based on new positions, determined by the location of the suspensions. If we consider these five units, the following methodology would be used: The first item must be the first failure; hence, it is assigned failure order number j = 1. The actual failure order number (or position) of the second failure, F₂, is in doubt. It could either be in position 2 or in position 3. Had S₁ not been withdrawn from the test at 9,500 hours, it could have operated successfully past 15,000 hours, thus placing F₂ in position 2. Alternatively, S₁ could also have failed before 15,000 hours, thus placing F₂ in position 3. In this case, the failure order number for F₂ will be some number between 2 and 3. To determine this number, consider the following:

We can find the number of ways the second failure can occur in either order number 2 (position 2) or order number 3 (position 3). The possible ways are listed next.

F2 in Position 2						OR	F2 in Position 3
1	2	3	4	5	6		1	2
F1	F1	F1	F1	F1	F1		F1	F1
F2	F2	F2	F2	F2	F2		S1	S1
S1	S2	F3	S1	S2	F3		F2	F2
S2	S1	S1	F3	F3	S2		S2	F3
F3	F3	S2	S2	S1	S1		F3	S2

It can be seen that F₂ can occur in the second position six ways and in the third position two ways. The most probable position is the average of these possible ways, or the mean order number ( MON ), given by:

[math]\displaystyle{ {{F}_{2}}=MO{{N}_{2}}=\frac{(6\times 2)+(2\times 3)}{6+2}=2.25 }[/math]

Using the same logic on the third failure, it can be located in position numbers 3, 4 and 5 in the possible ways listed next.

F3 in Position 3		OR	F3 in Position 4			OR	F3 in Position 5
1	2		1	2	3		1	2	3
F1	F1		F1	F1	F1		F1	F1	F1
F2	F2		S1	F2	F2		S1	F2	F2
F3	F3		F2	S1	S2		F2	S1	S2
S1	S2		F3	F3	F3		S2	S2	S1
S2	S1		S2	S2	S1		F3	F3	F3

Then, the mean order number for the third failure, (Item 5) is:

[math]\displaystyle{ MO{{N}_{3}}=\frac{(2\times 3)+(3\times 4)+(3\times 5)}{2+3+3}=4.125 }[/math]

Once the mean order number for each failure has been established, we obtain the median rank positions for these failures at their mean order number. Specifically, we obtain the median rank of the order numbers 1, 2.25 and 4.125 out of a sample size of 5, as given next.

Plotting Positions for the Failures(Sample Size=5)
Failure Number	MON	Median Rank Position(%)
1:F1	1	13%
2:F2	2.25	36%
3:F3	4.125	71%

Once the median rank values have been obtained, the probability plotting analysis is identical to that presented before. As you might have noticed, this methodology is rather laborious. Other techniques and shortcuts have been developed over the years to streamline this procedure. For more details on this method, see Kececioglu [20].

Shortfalls of this Rank Adjustment Method

Even though the rank adjustment method is the most widely used method for performing suspended items analysis, we would like to point out the following shortcoming. As you may have noticed from this analysis of suspended items, only the position where the failure occurred is taken into account, and not the exact time-to-suspension. For example, this methodology would yield the exact same results for the next two cases.

Case 1			Case 2
Item Number	State*"F" or "S"	Life of an item, hr	Item number	State*,"F" or "S"	Life of item, hr
1	F1	1,000	1	F1	1,000
2	S1	1,100	2	S1	9,700
3	S2	1,200	3	S2	9,800
4	S3	1,300	4	S3	9,900
5	F2	10,000	5	F2	10,000
* F − Failed, S − Suspended			* F − Failed, S − Suspended

This shortfall is significant when the number of failures is small and the number of suspensions is large and not spread uniformly between failures, as with these data. In cases like this, it is highly recommended that one use maximum likelihood estimation (MLE) to estimate the parameters instead of using least squares, since maximum likelihood does not look at ranks or plotting positions, but rather considers each unique time-to-failure or suspension. For the data given above, the results are as follows. The estimated parameters using the method just described are the same for both cases (1 and 2):

[math]\displaystyle{ \begin{array}{*{35}{l}} \widehat{\beta }= & \text{0}\text{.81} \\ \widehat{\eta }= & \text{11,417 hr} \\ \end{array} }[/math]

However, the MLE results for Case 1 are:

[math]\displaystyle{ \begin{array}{*{35}{l}} \widehat{\beta }= & \text{1}\text{.33} \\ \widehat{\eta }= & \text{6,900 hr} \\ \end{array} }[/math]

and the MLE results for Case 2 are:

[math]\displaystyle{ \begin{array}{*{35}{l}} \widehat{\beta }= & \text{0}\text{.9337} \\ \widehat{\eta }= & \text{21,348 hr} \\ \end{array} }[/math]

As we can see, there is a sizable difference in the results of the two sets calculated using MLE and the results using regression. The results for both cases are identical when using the regression estimation technique, as regression considers only the positions of the suspensions. The MLE results are quite different for the two cases, with the second case having a much larger value of η, which is due to the higher values of the suspension times in Case 2. This is because the maximum likelihood technique, unlike rank regression, considers the values of the suspensions when estimating the parameters. This is illustrated in the section MLE Analysis of Right Censored Data.

Probability Plots for Other Distributions

This same methodology can be applied to other distributions which have cdf equations that can be linearized. Different probability papers exist for each distribution, since different distributions have different cdf equations. Weibull++ automatically creates these plots for you when choosing a probability plot for a particular distribution. Special scales on these plots allow the parameter estimates to be derived directly from the plots, similar to the way β and η were obtained from the Weibull probability plot. These are discussed in the individual distribution chapters.

Let's assume six identical units are reliability tested at the same application and operation stress levels. All of these units fail during the test after operating for the following times (in hours): 96, 257, 498, 763, 1051 and 1744.

The steps for using the probability plotting method to determine the parameters of the exponential pdf representing the data are as follows:

Rank the times-to-failure in ascending order as shown next.

[math]\displaystyle{ \begin{matrix} \text{Failure} & \text{Failure Order Number} \\ \text{Time (Hr)} & \text{out of a Sample Size of 6} \\ \text{96} & \text{1} \\ \text{257} & \text{2} \\ \text{498} & \text{3} \\ \text{763} & \text{4} \\ \text{1,051} & \text{5} \\ \text{1,744} & \text{6} \\ \end{matrix}\,\! }[/math]

Obtain their median rank plotting positions. Median rank positions are used instead of other ranking methods because median ranks are at a specific confidence level (50%).

The times-to-failure, with their corresponding median ranks, are shown next:

[math]\displaystyle{ \begin{matrix} \text{Failure} & \text{Median} \\ \text{Time (Hr)} & \text{Rank, }% \\ \text{96} & \text{10}\text{.91} \\ \text{257} & \text{26}\text{.44} \\ \text{498} & \text{42}\text{.14} \\ \text{763} & \text{57}\text{.86} \\ \text{1,051} & \text{73}\text{.56} \\ \text{1,744} & \text{89}\text{.10} \\ \end{matrix}\,\! }[/math]

On an exponential probability paper, plot the times on the x-axis and their corresponding rank value on the y-axis. The next figure displays an example of an exponential probability paper. The paper is simply a log-linear paper.

Draw the best possible straight line that goes through the [math]\displaystyle{ t=0\,\! }[/math] and [math]\displaystyle{ (t)=100%\,\! }[/math] point and through the plotted points (as shown in the plot below).

At the [math]\displaystyle{ Q(t)=63.2%\,\! }[/math] or [math]\displaystyle{ R(t)=36.8%\,\! }[/math] ordinate point, draw a straight horizontal line until this line intersects the fitted straight line. Draw a vertical line through this intersection until it crosses the abscissa. The value at the intersection of the abscissa is the estimate of the mean. For this case, [math]\displaystyle{ \widehat{\mu }=833\,\! }[/math] hours which means that [math]\displaystyle{ \lambda =\tfrac{1}{\mu }=0.0012\,\! }[/math] (This is always at 63.2% because [math]\displaystyle{ (T)=1-{{e}^{-\tfrac{\mu }{\mu }}}=1-{{e}^{-1}}=0.632=63.2%)\,\! }[/math].

Now any reliability value for any mission time [math]\displaystyle{ t\,\! }[/math] can be obtained. For example, the reliability for a mission of 15 hours, or any other time, can now be obtained either from the plot or analytically.

To obtain the value from the plot, draw a vertical line from the abscissa, at [math]\displaystyle{ t=15\,\! }[/math] hours, to the fitted line. Draw a horizontal line from this intersection to the ordinate and read [math]\displaystyle{ R(t)\,\! }[/math]. In this case, [math]\displaystyle{ R(t=15)=98.15%\,\! }[/math]. This can also be obtained analytically, from the exponential reliability function.

Some Shortfalls of Manual Probability Plotting

Besides the most obvious drawback to probability plotting, which is the amount of effort required, manual probability plotting is not always consistent in the results. Two people plotting a straight line through a set of points will not always draw this line the same way, and thus will come up with slightly different results. This method was used primarily before the widespread use of computers that could easily perform the calculations for more complicated parameter estimation methods, such as the least squares and maximum likelihood methods.

Parameter estimation refers to the process of using sample data (in our case times-to-failure or success data) to estimate the parameters of the selected distribution. Several parameter estimation methods are available. This section presents an overview of the available methods used in life data analysis. More specifically, we start with the relatively simple method of probability plotting and continue with the more sophisticated methods of Rank Regression (or Least Squares) and Maximum Likelihood.

Probability Plotting

The least mathematically intensive method for parameter estimation is the method of probability plotting. As the term implies, probability plotting involves a physical plot of the data on specially constructed probability plotting paper. This method is easily implemented by hand, given that one can obtain the appropriate probability plotting paper.

Illustrating the Method for the 2-Parameter Weibull Distribution

The method of probability plotting takes the cdf of the distribution and attempts to linearize it by employing a specially constructed paper. This is best illustrated using the 2-parameter Weibull distribution.

In the case of the two-parameter Weibull distribution, the cdf (also the unreliability Q(t)) is given by:

[math]\displaystyle{ F(t)=Q(t)=1-{e^{-\left(\tfrac{t}{\eta}\right)^{\beta}}} }[/math]

Linearizing the Weibull Unreliability Function
This function can then be linearized (i.e., put in the common form of y = m'x + b format) as follows:

[math]\displaystyle{ \begin{align} Q(t)= & 1-{e^{-\left(\tfrac{t}{\eta}\right)^{\beta}}} \\ \ln (1-Q(t))= & \ln \left[ {e^{-\left(\tfrac{t}{\eta}\right)^{\beta}}} \right] \\ \ln (1-Q(t))=& -\left(\tfrac{t}{\eta}\right)^{\beta} \\ \ln ( -\ln (1-Q(t)))= & \beta \left(\ln \left( \frac{t}{\eta }\right)\right) \\ \ln \left( \ln \left( \frac{1}{1-Q(t)}\right) \right) = & \beta\ln{ t} -\beta(\eta ) \\ \end{align} }[/math]

Then by setting

[math]\displaystyle{ y=\ln \left( \ln \left( \frac{1}{1-Q(t)} \right) \right) }[/math]

and

[math]\displaystyle{ x=\ln \left( t \right) }[/math]

the equation can then be rewritten as:

[math]\displaystyle{ y=\beta x-\beta \ln \left( \eta \right) }[/math]

which is now a linear equation with a slope of

Slope = m = β

and an intercept of

Intercept[math]\displaystyle{ =b=-\beta \cdot ln(\eta) }[/math].

Constructing the Paper

The next task is to construct the Weibull probability plotting paper with the appropriate y - and x-axes. The x-axis trasnformation is simply logarithmic. The y-axis, is a bit more complex requiring a double log reciprocal transformation, or 

[math]\displaystyle{ y=\ln \left(\ln \left( \frac{1}{1-Q(t)} ) \right) \right) }[/math]

where Q(t) is the unreliability.

Such papers have been created by different vendors and are called probability plotting papers. Weibull.com has different plotting papers available for download.

To illustrate, consider the following probability plot on a slightly different type of Weibull probability paper.

This paper is constructed based on the mentioned y - and x - transformations, where the y-axis represents unreliability and the x-axis represents time. Both of these values must be known for each time-to-failure point we want to plot.

Then, given the y and x value for each point, the points can easily be put on the plot. Once the points have been placed on the plot, the best possible straight line is drawn through these points. Once the line has been drawn, the slope of the line can be obtained (some probability papers include a slope indicator to simplify this calculation). This is the parameter β, which is the value of the slope. To determine the scale parameter, η (also called the characteristic life), one reads the time from the x-axis corresponding to Q(t)=63.2%.

Note that at

[math]\displaystyle{ \begin{align} Q(t=\eta)= & 1-{{e}^{-{{\left( \tfrac{t}{\eta } \right)}^{\beta }}}} \\ = & 1-{{e}^{-1}} \\ = & 0.632 \\ = & 63.2% \end{align} }[/math]

Thus, if we enter the y axis at Q(t)=63.2%, the corresponding value of t will be equal to η. Thus, using this simple methodology, the parameters of the Weibull distribution can be estimated.

Determining the x and y Position of the Plot Points

The points on the plot represent our data or, more specifically, our times-to-failure data. If, for example, we tested four units that failed at 10, 20, 30 and 40 hours, then we would use these times as our x values or time values.

Determining the appropriate y plotting positions, or the unreliability values, is a little more complex. To determine the y plotting positions, we must first determine a value indicating the corresponding unreliability for that failure. In other words, we need to obtain the cumulative percent failed for each time-to-failure. For example, the cumulative percent failed by 10 hours may be 25%, by 20 hours 50%, and so forth. This is a simple method illustrating the idea. The problem with this simple method is the fact that the 100% point is not defined on most probability plots; thus, an alternative and more robust approach must be used. The most widely used method of determining this value is the method of obtaining the median rank for each failure. This is discussed next.

Median Ranks

Median ranks are used to obtain an estimate of the unreliability for each failure. It is the value that the true probability of failure, Q(T_j), should have at the j^t'h failure out of a sample of N units at the 50% confidence level.

The rank can be found for any percentage point, P, greater than zero and less than one, by solving the cumulative binomial equation for Z . This represents the rank, or unreliability estimate, for the j^t'h failure in the following equation for the cumulative binomial:

[math]\displaystyle{ P=\underset{k=j}{\overset{N}{\mathop \sum }}\,\left( \begin{matrix} N \\ k \\ \end{matrix} \right){{Z}^{k}}{{\left( 1-Z \right)}^{N-k}} }[/math]

where N is the sample size and j the order number.

The median rank is obtained by solving this equation for Z at P = 0.50,

[math]\displaystyle{ 0.50=\underset{k=j}{\overset{N}{\mathop \sum }}\,\left( \begin{matrix} N \\ k \\ \end{matrix} \right){{Z}^{k}}{{\left( 1-Z \right)}^{N-k}} }[/math]

For example, if N=4 and we have four failures, we would solve the median rank equation for the value of Z four times; once for each failure with j= 1, 2, 3 and 4. This result can then be used as the unreliability estimate for each failure or the y plotting position. (See also the chapter The Weibull distribution for a step-by-step example of this method.) The solution of cumuative binomial equation for Z requires the use of numerical methods.

Beta and F Distributions Approach

A more straightforward and easier method of estimating median ranks is by applying two transformations to the cumulative binomial equation, first to the beta distribution and then to the F distribution, resulting in [12, 13],

[math]\displaystyle{ \begin{array}{*{35}{l}} MR & = & \tfrac{1}{1+\tfrac{N-j+1}{j}{{F}_{0.50;m;n}}} \\ m & = & 2(N-j+1) \\ n & = & 2j \\ \end{array} }[/math]

where F_0.50;m;n denotes the F distribution at the 0.50 point, with m and n degrees of freedom, for failure j out of N units.

Benard's Approximation for Median Ranks

Another quick, and less accurate, approximation of the median ranks is also given by:

[math]\displaystyle{ MR = \frac{{j - 0.3}}{{N + 0.4}} }[/math]

This approximation of the median ranks is also known as Benard's approximation.

Kaplan-Meier

The Kaplan-Meier estimator (also known as the product limit estimator) is used as an alternative to the median ranks method for calculating the estimates of the unreliability for probability plotting purposes. The equation of the estimator is given by:

[math]\displaystyle{ \widehat{F}({{t}_{i}})=1-\underset{j=1}{\overset{i}{\mathop \prod }}\,\frac{{{n}_{j}}-{{r}_{j}}}{{{n}_{j}}},\text{ }i=1,...,m }[/math]

where:

[math]\displaystyle{ \begin{align} m = & {\text{total number of data points}} \\ n = & {\text{the total number of units}} \\ {n_i} = & n - \sum_{j = 0}^{i - 1}{s_j} - \sum_{j = 0}^{i - 1}{r_j}, \text{i = 1,...,m }\\ {r_j} = & {\text{ number of failures in the }}{j^{th}}{\text{ data group, and}} \\ {s_j} = & {\text{number of surviving units in the }}{j^{th}}{\text{ data group}} \\ \end{align} }[/math]

Analysis of Right Censored (Suspended) Data

All available data should be considered in the analysis of times-to-failure data. This includes the case when a particular unit in a sample has been removed from the test prior to failure. An item, or unit, which is removed from a reliability test prior to failure, or a unit which is in the field and is still operating at the time the reliability of these units is to be determined, is called a suspended item or right censored observation or right censored data point. Suspended items analysis would also be considered when:

1. We need to make an analysis of the available results before test completion.
2. The failure modes which are occurring are different than those anticipated and such units are withdrawn from the test.
3. We need to analyze a single mode and the actual data set comprises multiple modes.
4. A warranty analysis is to be made of all units in the field (non-failed and failed units). The non-failed units are considered to be suspended items (or right censored).

Probability Plotting and Rank Regression Analysis of Right Censored Data

When using the probability plotting or rank regression method to accommodate the fact that units in the data set did not fail, or were suspended, we need to adjust their probability of failure, or unreliability. As discussed before, estimates of the unreliability for complete data are obtained using the median ranks approach. The following methodology illustrates how adjusted median ranks are computed to account for right censored data. To better illustrate the methodology, consider the following example in  [20] where five items are tested resulting in three failures and two suspensions.

Item Number(Position)	State*"F" or "S"	Life of item, hr
1	F1	5,100
2	S1	9,500
3	F2	15,000
4	S2	22,000
5	F3	40,000

The methodology for plotting suspended items involves adjusting the rank positions and plotting the data based on new positions, determined by the location of the suspensions. If we consider these five units, the following methodology would be used: The first item must be the first failure; hence, it is assigned failure order number j = 1. The actual failure order number (or position) of the second failure, F₂, is in doubt. It could either be in position 2 or in position 3. Had S₁ not been withdrawn from the test at 9,500 hours, it could have operated successfully past 15,000 hours, thus placing F₂ in position 2. Alternatively, S₁ could also have failed before 15,000 hours, thus placing F₂ in position 3. In this case, the failure order number for F₂ will be some number between 2 and 3. To determine this number, consider the following:

We can find the number of ways the second failure can occur in either order number 2 (position 2) or order number 3 (position 3). The possible ways are listed next.

F2 in Position 2						OR	F2 in Position 3
1	2	3	4	5	6		1	2
F1	F1	F1	F1	F1	F1		F1	F1
F2	F2	F2	F2	F2	F2		S1	S1
S1	S2	F3	S1	S2	F3		F2	F2
S2	S1	S1	F3	F3	S2		S2	F3
F3	F3	S2	S2	S1	S1		F3	S2

It can be seen that F₂ can occur in the second position six ways and in the third position two ways. The most probable position is the average of these possible ways, or the mean order number ( MON ), given by:

[math]\displaystyle{ {{F}_{2}}=MO{{N}_{2}}=\frac{(6\times 2)+(2\times 3)}{6+2}=2.25 }[/math]

Using the same logic on the third failure, it can be located in position numbers 3, 4 and 5 in the possible ways listed next.

F3 in Position 3		OR	F3 in Position 4			OR	F3 in Position 5
1	2		1	2	3		1	2	3
F1	F1		F1	F1	F1		F1	F1	F1
F2	F2		S1	F2	F2		S1	F2	F2
F3	F3		F2	S1	S2		F2	S1	S2
S1	S2		F3	F3	F3		S2	S2	S1
S2	S1		S2	S2	S1		F3	F3	F3

Then, the mean order number for the third failure, (Item 5) is:

[math]\displaystyle{ MO{{N}_{3}}=\frac{(2\times 3)+(3\times 4)+(3\times 5)}{2+3+3}=4.125 }[/math]

Once the mean order number for each failure has been established, we obtain the median rank positions for these failures at their mean order number. Specifically, we obtain the median rank of the order numbers 1, 2.25 and 4.125 out of a sample size of 5, as given next.

Plotting Positions for the Failures(Sample Size=5)
Failure Number	MON	Median Rank Position(%)
1:F1	1	13%
2:F2	2.25	36%
3:F3	4.125	71%

Once the median rank values have been obtained, the probability plotting analysis is identical to that presented before. As you might have noticed, this methodology is rather laborious. Other techniques and shortcuts have been developed over the years to streamline this procedure. For more details on this method, see Kececioglu [20].

Shortfalls of this Rank Adjustment Method

Even though the rank adjustment method is the most widely used method for performing suspended items analysis, we would like to point out the following shortcoming. As you may have noticed from this analysis of suspended items, only the position where the failure occurred is taken into account, and not the exact time-to-suspension. For example, this methodology would yield the exact same results for the next two cases.

Case 1			Case 2
Item Number	State*"F" or "S"	Life of an item, hr	Item number	State*,"F" or "S"	Life of item, hr
1	F1	1,000	1	F1	1,000
2	S1	1,100	2	S1	9,700
3	S2	1,200	3	S2	9,800
4	S3	1,300	4	S3	9,900
5	F2	10,000	5	F2	10,000
* F − Failed, S − Suspended			* F − Failed, S − Suspended

This shortfall is significant when the number of failures is small and the number of suspensions is large and not spread uniformly between failures, as with these data. In cases like this, it is highly recommended that one use maximum likelihood estimation (MLE) to estimate the parameters instead of using least squares, since maximum likelihood does not look at ranks or plotting positions, but rather considers each unique time-to-failure or suspension. For the data given above, the results are as follows. The estimated parameters using the method just described are the same for both cases (1 and 2):

[math]\displaystyle{ \begin{array}{*{35}{l}} \widehat{\beta }= & \text{0}\text{.81} \\ \widehat{\eta }= & \text{11,417 hr} \\ \end{array} }[/math]

However, the MLE results for Case 1 are:

[math]\displaystyle{ \begin{array}{*{35}{l}} \widehat{\beta }= & \text{1}\text{.33} \\ \widehat{\eta }= & \text{6,900 hr} \\ \end{array} }[/math]

and the MLE results for Case 2 are:

[math]\displaystyle{ \begin{array}{*{35}{l}} \widehat{\beta }= & \text{0}\text{.9337} \\ \widehat{\eta }= & \text{21,348 hr} \\ \end{array} }[/math]

As we can see, there is a sizable difference in the results of the two sets calculated using MLE and the results using regression. The results for both cases are identical when using the regression estimation technique, as regression considers only the positions of the suspensions. The MLE results are quite different for the two cases, with the second case having a much larger value of η, which is due to the higher values of the suspension times in Case 2. This is because the maximum likelihood technique, unlike rank regression, considers the values of the suspensions when estimating the parameters. This is illustrated in the section MLE Analysis of Right Censored Data.

Probability Plots for Other Distributions

This same methodology can be applied to other distributions which have cdf equations that can be linearized. Different probability papers exist for each distribution, since different distributions have different cdf equations. Weibull++ automatically creates these plots for you when choosing a probability plot for a particular distribution. Special scales on these plots allow the parameter estimates to be derived directly from the plots, similar to the way β and η were obtained from the Weibull probability plot. These are discussed in the individual distribution chapters.

Let's assume six identical units are reliability tested at the same application and operation stress levels. All of these units fail during the test after operating for the following times (in hours): 96, 257, 498, 763, 1051 and 1744.

The steps for using the probability plotting method to determine the parameters of the exponential pdf representing the data are as follows:

Rank the times-to-failure in ascending order as shown next.

[math]\displaystyle{ \begin{matrix} \text{Failure} & \text{Failure Order Number} \\ \text{Time (Hr)} & \text{out of a Sample Size of 6} \\ \text{96} & \text{1} \\ \text{257} & \text{2} \\ \text{498} & \text{3} \\ \text{763} & \text{4} \\ \text{1,051} & \text{5} \\ \text{1,744} & \text{6} \\ \end{matrix}\,\! }[/math]

Obtain their median rank plotting positions. Median rank positions are used instead of other ranking methods because median ranks are at a specific confidence level (50%).

The times-to-failure, with their corresponding median ranks, are shown next:

[math]\displaystyle{ \begin{matrix} \text{Failure} & \text{Median} \\ \text{Time (Hr)} & \text{Rank, }% \\ \text{96} & \text{10}\text{.91} \\ \text{257} & \text{26}\text{.44} \\ \text{498} & \text{42}\text{.14} \\ \text{763} & \text{57}\text{.86} \\ \text{1,051} & \text{73}\text{.56} \\ \text{1,744} & \text{89}\text{.10} \\ \end{matrix}\,\! }[/math]

On an exponential probability paper, plot the times on the x-axis and their corresponding rank value on the y-axis. The next figure displays an example of an exponential probability paper. The paper is simply a log-linear paper.

Draw the best possible straight line that goes through the [math]\displaystyle{ t=0\,\! }[/math] and [math]\displaystyle{ (t)=100%\,\! }[/math] point and through the plotted points (as shown in the plot below).

At the [math]\displaystyle{ Q(t)=63.2%\,\! }[/math] or [math]\displaystyle{ R(t)=36.8%\,\! }[/math] ordinate point, draw a straight horizontal line until this line intersects the fitted straight line. Draw a vertical line through this intersection until it crosses the abscissa. The value at the intersection of the abscissa is the estimate of the mean. For this case, [math]\displaystyle{ \widehat{\mu }=833\,\! }[/math] hours which means that [math]\displaystyle{ \lambda =\tfrac{1}{\mu }=0.0012\,\! }[/math] (This is always at 63.2% because [math]\displaystyle{ (T)=1-{{e}^{-\tfrac{\mu }{\mu }}}=1-{{e}^{-1}}=0.632=63.2%)\,\! }[/math].

Now any reliability value for any mission time [math]\displaystyle{ t\,\! }[/math] can be obtained. For example, the reliability for a mission of 15 hours, or any other time, can now be obtained either from the plot or analytically.

To obtain the value from the plot, draw a vertical line from the abscissa, at [math]\displaystyle{ t=15\,\! }[/math] hours, to the fitted line. Draw a horizontal line from this intersection to the ordinate and read [math]\displaystyle{ R(t)\,\! }[/math]. In this case, [math]\displaystyle{ R(t=15)=98.15%\,\! }[/math]. This can also be obtained analytically, from the exponential reliability function.

Some Shortfalls of Manual Probability Plotting

Besides the most obvious drawback to probability plotting, which is the amount of effort required, manual probability plotting is not always consistent in the results. Two people plotting a straight line through a set of points will not always draw this line the same way, and thus will come up with slightly different results. This method was used primarily before the widespread use of computers that could easily perform the calculations for more complicated parameter estimation methods, such as the least squares and maximum likelihood methods.

Parameter estimation refers to the process of using sample data (in our case times-to-failure or success data) to estimate the parameters of the selected distribution. Several parameter estimation methods are available. This section presents an overview of the available methods used in life data analysis. More specifically, we start with the relatively simple method of probability plotting and continue with the more sophisticated methods of Rank Regression (or Least Squares) and Maximum Likelihood.

Probability Plotting

The least mathematically intensive method for parameter estimation is the method of probability plotting. As the term implies, probability plotting involves a physical plot of the data on specially constructed probability plotting paper. This method is easily implemented by hand, given that one can obtain the appropriate probability plotting paper.

Illustrating the Method for the 2-Parameter Weibull Distribution

The method of probability plotting takes the cdf of the distribution and attempts to linearize it by employing a specially constructed paper. This is best illustrated using the 2-parameter Weibull distribution.

In the case of the two-parameter Weibull distribution, the cdf (also the unreliability Q(t)) is given by:

[math]\displaystyle{ F(t)=Q(t)=1-{e^{-\left(\tfrac{t}{\eta}\right)^{\beta}}} }[/math]

Linearizing the Weibull Unreliability Function
This function can then be linearized (i.e., put in the common form of y = m'x + b format) as follows:

[math]\displaystyle{ \begin{align} Q(t)= & 1-{e^{-\left(\tfrac{t}{\eta}\right)^{\beta}}} \\ \ln (1-Q(t))= & \ln \left[ {e^{-\left(\tfrac{t}{\eta}\right)^{\beta}}} \right] \\ \ln (1-Q(t))=& -\left(\tfrac{t}{\eta}\right)^{\beta} \\ \ln ( -\ln (1-Q(t)))= & \beta \left(\ln \left( \frac{t}{\eta }\right)\right) \\ \ln \left( \ln \left( \frac{1}{1-Q(t)}\right) \right) = & \beta\ln{ t} -\beta(\eta ) \\ \end{align} }[/math]

Then by setting

[math]\displaystyle{ y=\ln \left( \ln \left( \frac{1}{1-Q(t)} \right) \right) }[/math]

and

[math]\displaystyle{ x=\ln \left( t \right) }[/math]

the equation can then be rewritten as:

[math]\displaystyle{ y=\beta x-\beta \ln \left( \eta \right) }[/math]

which is now a linear equation with a slope of

Slope = m = β

and an intercept of

Intercept[math]\displaystyle{ =b=-\beta \cdot ln(\eta) }[/math].

Constructing the Paper

The next task is to construct the Weibull probability plotting paper with the appropriate y - and x-axes. The x-axis trasnformation is simply logarithmic. The y-axis, is a bit more complex requiring a double log reciprocal transformation, or 

[math]\displaystyle{ y=\ln \left(\ln \left( \frac{1}{1-Q(t)} ) \right) \right) }[/math]

where Q(t) is the unreliability.

Such papers have been created by different vendors and are called probability plotting papers. Weibull.com has different plotting papers available for download.

To illustrate, consider the following probability plot on a slightly different type of Weibull probability paper.

This paper is constructed based on the mentioned y - and x - transformations, where the y-axis represents unreliability and the x-axis represents time. Both of these values must be known for each time-to-failure point we want to plot.

Then, given the y and x value for each point, the points can easily be put on the plot. Once the points have been placed on the plot, the best possible straight line is drawn through these points. Once the line has been drawn, the slope of the line can be obtained (some probability papers include a slope indicator to simplify this calculation). This is the parameter β, which is the value of the slope. To determine the scale parameter, η (also called the characteristic life), one reads the time from the x-axis corresponding to Q(t)=63.2%.

Note that at

[math]\displaystyle{ \begin{align} Q(t=\eta)= & 1-{{e}^{-{{\left( \tfrac{t}{\eta } \right)}^{\beta }}}} \\ = & 1-{{e}^{-1}} \\ = & 0.632 \\ = & 63.2% \end{align} }[/math]

Thus, if we enter the y axis at Q(t)=63.2%, the corresponding value of t will be equal to η. Thus, using this simple methodology, the parameters of the Weibull distribution can be estimated.

Determining the x and y Position of the Plot Points

The points on the plot represent our data or, more specifically, our times-to-failure data. If, for example, we tested four units that failed at 10, 20, 30 and 40 hours, then we would use these times as our x values or time values.

Determining the appropriate y plotting positions, or the unreliability values, is a little more complex. To determine the y plotting positions, we must first determine a value indicating the corresponding unreliability for that failure. In other words, we need to obtain the cumulative percent failed for each time-to-failure. For example, the cumulative percent failed by 10 hours may be 25%, by 20 hours 50%, and so forth. This is a simple method illustrating the idea. The problem with this simple method is the fact that the 100% point is not defined on most probability plots; thus, an alternative and more robust approach must be used. The most widely used method of determining this value is the method of obtaining the median rank for each failure. This is discussed next.

Median Ranks

Median ranks are used to obtain an estimate of the unreliability for each failure. It is the value that the true probability of failure, Q(T_j), should have at the j^t'h failure out of a sample of N units at the 50% confidence level.

The rank can be found for any percentage point, P, greater than zero and less than one, by solving the cumulative binomial equation for Z . This represents the rank, or unreliability estimate, for the j^t'h failure in the following equation for the cumulative binomial:

[math]\displaystyle{ P=\underset{k=j}{\overset{N}{\mathop \sum }}\,\left( \begin{matrix} N \\ k \\ \end{matrix} \right){{Z}^{k}}{{\left( 1-Z \right)}^{N-k}} }[/math]

where N is the sample size and j the order number.

The median rank is obtained by solving this equation for Z at P = 0.50,

[math]\displaystyle{ 0.50=\underset{k=j}{\overset{N}{\mathop \sum }}\,\left( \begin{matrix} N \\ k \\ \end{matrix} \right){{Z}^{k}}{{\left( 1-Z \right)}^{N-k}} }[/math]

For example, if N=4 and we have four failures, we would solve the median rank equation for the value of Z four times; once for each failure with j= 1, 2, 3 and 4. This result can then be used as the unreliability estimate for each failure or the y plotting position. (See also the chapter The Weibull distribution for a step-by-step example of this method.) The solution of cumuative binomial equation for Z requires the use of numerical methods.

Beta and F Distributions Approach

A more straightforward and easier method of estimating median ranks is by applying two transformations to the cumulative binomial equation, first to the beta distribution and then to the F distribution, resulting in [12, 13],

[math]\displaystyle{ \begin{array}{*{35}{l}} MR & = & \tfrac{1}{1+\tfrac{N-j+1}{j}{{F}_{0.50;m;n}}} \\ m & = & 2(N-j+1) \\ n & = & 2j \\ \end{array} }[/math]

where F_0.50;m;n denotes the F distribution at the 0.50 point, with m and n degrees of freedom, for failure j out of N units.

Benard's Approximation for Median Ranks

Another quick, and less accurate, approximation of the median ranks is also given by:

[math]\displaystyle{ MR = \frac{{j - 0.3}}{{N + 0.4}} }[/math]

This approximation of the median ranks is also known as Benard's approximation.

Kaplan-Meier

The Kaplan-Meier estimator (also known as the product limit estimator) is used as an alternative to the median ranks method for calculating the estimates of the unreliability for probability plotting purposes. The equation of the estimator is given by:

[math]\displaystyle{ \widehat{F}({{t}_{i}})=1-\underset{j=1}{\overset{i}{\mathop \prod }}\,\frac{{{n}_{j}}-{{r}_{j}}}{{{n}_{j}}},\text{ }i=1,...,m }[/math]

where:

[math]\displaystyle{ \begin{align} m = & {\text{total number of data points}} \\ n = & {\text{the total number of units}} \\ {n_i} = & n - \sum_{j = 0}^{i - 1}{s_j} - \sum_{j = 0}^{i - 1}{r_j}, \text{i = 1,...,m }\\ {r_j} = & {\text{ number of failures in the }}{j^{th}}{\text{ data group, and}} \\ {s_j} = & {\text{number of surviving units in the }}{j^{th}}{\text{ data group}} \\ \end{align} }[/math]

Analysis of Right Censored (Suspended) Data

All available data should be considered in the analysis of times-to-failure data. This includes the case when a particular unit in a sample has been removed from the test prior to failure. An item, or unit, which is removed from a reliability test prior to failure, or a unit which is in the field and is still operating at the time the reliability of these units is to be determined, is called a suspended item or right censored observation or right censored data point. Suspended items analysis would also be considered when:

1. We need to make an analysis of the available results before test completion.
2. The failure modes which are occurring are different than those anticipated and such units are withdrawn from the test.
3. We need to analyze a single mode and the actual data set comprises multiple modes.
4. A warranty analysis is to be made of all units in the field (non-failed and failed units). The non-failed units are considered to be suspended items (or right censored).

Probability Plotting and Rank Regression Analysis of Right Censored Data

When using the probability plotting or rank regression method to accommodate the fact that units in the data set did not fail, or were suspended, we need to adjust their probability of failure, or unreliability. As discussed before, estimates of the unreliability for complete data are obtained using the median ranks approach. The following methodology illustrates how adjusted median ranks are computed to account for right censored data. To better illustrate the methodology, consider the following example in  [20] where five items are tested resulting in three failures and two suspensions.

Item Number(Position)	State*"F" or "S"	Life of item, hr
1	F1	5,100
2	S1	9,500
3	F2	15,000
4	S2	22,000
5	F3	40,000

The methodology for plotting suspended items involves adjusting the rank positions and plotting the data based on new positions, determined by the location of the suspensions. If we consider these five units, the following methodology would be used: The first item must be the first failure; hence, it is assigned failure order number j = 1. The actual failure order number (or position) of the second failure, F₂, is in doubt. It could either be in position 2 or in position 3. Had S₁ not been withdrawn from the test at 9,500 hours, it could have operated successfully past 15,000 hours, thus placing F₂ in position 2. Alternatively, S₁ could also have failed before 15,000 hours, thus placing F₂ in position 3. In this case, the failure order number for F₂ will be some number between 2 and 3. To determine this number, consider the following:

We can find the number of ways the second failure can occur in either order number 2 (position 2) or order number 3 (position 3). The possible ways are listed next.

F2 in Position 2						OR	F2 in Position 3
1	2	3	4	5	6		1	2
F1	F1	F1	F1	F1	F1		F1	F1
F2	F2	F2	F2	F2	F2		S1	S1
S1	S2	F3	S1	S2	F3		F2	F2
S2	S1	S1	F3	F3	S2		S2	F3
F3	F3	S2	S2	S1	S1		F3	S2

It can be seen that F₂ can occur in the second position six ways and in the third position two ways. The most probable position is the average of these possible ways, or the mean order number ( MON ), given by:

[math]\displaystyle{ {{F}_{2}}=MO{{N}_{2}}=\frac{(6\times 2)+(2\times 3)}{6+2}=2.25 }[/math]

Using the same logic on the third failure, it can be located in position numbers 3, 4 and 5 in the possible ways listed next.

F3 in Position 3		OR	F3 in Position 4			OR	F3 in Position 5
1	2		1	2	3		1	2	3
F1	F1		F1	F1	F1		F1	F1	F1
F2	F2		S1	F2	F2		S1	F2	F2
F3	F3		F2	S1	S2		F2	S1	S2
S1	S2		F3	F3	F3		S2	S2	S1
S2	S1		S2	S2	S1		F3	F3	F3

Then, the mean order number for the third failure, (Item 5) is:

[math]\displaystyle{ MO{{N}_{3}}=\frac{(2\times 3)+(3\times 4)+(3\times 5)}{2+3+3}=4.125 }[/math]

Once the mean order number for each failure has been established, we obtain the median rank positions for these failures at their mean order number. Specifically, we obtain the median rank of the order numbers 1, 2.25 and 4.125 out of a sample size of 5, as given next.

Plotting Positions for the Failures(Sample Size=5)
Failure Number	MON	Median Rank Position(%)
1:F1	1	13%
2:F2	2.25	36%
3:F3	4.125	71%

Once the median rank values have been obtained, the probability plotting analysis is identical to that presented before. As you might have noticed, this methodology is rather laborious. Other techniques and shortcuts have been developed over the years to streamline this procedure. For more details on this method, see Kececioglu [20].

Shortfalls of this Rank Adjustment Method

Even though the rank adjustment method is the most widely used method for performing suspended items analysis, we would like to point out the following shortcoming. As you may have noticed from this analysis of suspended items, only the position where the failure occurred is taken into account, and not the exact time-to-suspension. For example, this methodology would yield the exact same results for the next two cases.

Case 1			Case 2
Item Number	State*"F" or "S"	Life of an item, hr	Item number	State*,"F" or "S"	Life of item, hr
1	F1	1,000	1	F1	1,000
2	S1	1,100	2	S1	9,700
3	S2	1,200	3	S2	9,800
4	S3	1,300	4	S3	9,900
5	F2	10,000	5	F2	10,000
* F − Failed, S − Suspended			* F − Failed, S − Suspended

This shortfall is significant when the number of failures is small and the number of suspensions is large and not spread uniformly between failures, as with these data. In cases like this, it is highly recommended that one use maximum likelihood estimation (MLE) to estimate the parameters instead of using least squares, since maximum likelihood does not look at ranks or plotting positions, but rather considers each unique time-to-failure or suspension. For the data given above, the results are as follows. The estimated parameters using the method just described are the same for both cases (1 and 2):

[math]\displaystyle{ \begin{array}{*{35}{l}} \widehat{\beta }= & \text{0}\text{.81} \\ \widehat{\eta }= & \text{11,417 hr} \\ \end{array} }[/math]

However, the MLE results for Case 1 are:

[math]\displaystyle{ \begin{array}{*{35}{l}} \widehat{\beta }= & \text{1}\text{.33} \\ \widehat{\eta }= & \text{6,900 hr} \\ \end{array} }[/math]

and the MLE results for Case 2 are:

[math]\displaystyle{ \begin{array}{*{35}{l}} \widehat{\beta }= & \text{0}\text{.9337} \\ \widehat{\eta }= & \text{21,348 hr} \\ \end{array} }[/math]

As we can see, there is a sizable difference in the results of the two sets calculated using MLE and the results using regression. The results for both cases are identical when using the regression estimation technique, as regression considers only the positions of the suspensions. The MLE results are quite different for the two cases, with the second case having a much larger value of η, which is due to the higher values of the suspension times in Case 2. This is because the maximum likelihood technique, unlike rank regression, considers the values of the suspensions when estimating the parameters. This is illustrated in the section MLE Analysis of Right Censored Data.

Probability Plots for Other Distributions

This same methodology can be applied to other distributions which have cdf equations that can be linearized. Different probability papers exist for each distribution, since different distributions have different cdf equations. Weibull++ automatically creates these plots for you when choosing a probability plot for a particular distribution. Special scales on these plots allow the parameter estimates to be derived directly from the plots, similar to the way β and η were obtained from the Weibull probability plot. These are discussed in the individual distribution chapters.

Let's assume six identical units are reliability tested at the same application and operation stress levels. All of these units fail during the test after operating for the following times (in hours): 96, 257, 498, 763, 1051 and 1744.

The steps for using the probability plotting method to determine the parameters of the exponential pdf representing the data are as follows:

Rank the times-to-failure in ascending order as shown next.

[math]\displaystyle{ \begin{matrix} \text{Failure} & \text{Failure Order Number} \\ \text{Time (Hr)} & \text{out of a Sample Size of 6} \\ \text{96} & \text{1} \\ \text{257} & \text{2} \\ \text{498} & \text{3} \\ \text{763} & \text{4} \\ \text{1,051} & \text{5} \\ \text{1,744} & \text{6} \\ \end{matrix}\,\! }[/math]

Obtain their median rank plotting positions. Median rank positions are used instead of other ranking methods because median ranks are at a specific confidence level (50%).

The times-to-failure, with their corresponding median ranks, are shown next:

[math]\displaystyle{ \begin{matrix} \text{Failure} & \text{Median} \\ \text{Time (Hr)} & \text{Rank, }% \\ \text{96} & \text{10}\text{.91} \\ \text{257} & \text{26}\text{.44} \\ \text{498} & \text{42}\text{.14} \\ \text{763} & \text{57}\text{.86} \\ \text{1,051} & \text{73}\text{.56} \\ \text{1,744} & \text{89}\text{.10} \\ \end{matrix}\,\! }[/math]

On an exponential probability paper, plot the times on the x-axis and their corresponding rank value on the y-axis. The next figure displays an example of an exponential probability paper. The paper is simply a log-linear paper.

Draw the best possible straight line that goes through the [math]\displaystyle{ t=0\,\! }[/math] and [math]\displaystyle{ (t)=100%\,\! }[/math] point and through the plotted points (as shown in the plot below).

At the [math]\displaystyle{ Q(t)=63.2%\,\! }[/math] or [math]\displaystyle{ R(t)=36.8%\,\! }[/math] ordinate point, draw a straight horizontal line until this line intersects the fitted straight line. Draw a vertical line through this intersection until it crosses the abscissa. The value at the intersection of the abscissa is the estimate of the mean. For this case, [math]\displaystyle{ \widehat{\mu }=833\,\! }[/math] hours which means that [math]\displaystyle{ \lambda =\tfrac{1}{\mu }=0.0012\,\! }[/math] (This is always at 63.2% because [math]\displaystyle{ (T)=1-{{e}^{-\tfrac{\mu }{\mu }}}=1-{{e}^{-1}}=0.632=63.2%)\,\! }[/math].

Now any reliability value for any mission time [math]\displaystyle{ t\,\! }[/math] can be obtained. For example, the reliability for a mission of 15 hours, or any other time, can now be obtained either from the plot or analytically.

To obtain the value from the plot, draw a vertical line from the abscissa, at [math]\displaystyle{ t=15\,\! }[/math] hours, to the fitted line. Draw a horizontal line from this intersection to the ordinate and read [math]\displaystyle{ R(t)\,\! }[/math]. In this case, [math]\displaystyle{ R(t=15)=98.15%\,\! }[/math]. This can also be obtained analytically, from the exponential reliability function.

Some Shortfalls of Manual Probability Plotting

Besides the most obvious drawback to probability plotting, which is the amount of effort required, manual probability plotting is not always consistent in the results. Two people plotting a straight line through a set of points will not always draw this line the same way, and thus will come up with slightly different results. This method was used primarily before the widespread use of computers that could easily perform the calculations for more complicated parameter estimation methods, such as the least squares and maximum likelihood methods.