Reliability Data Analysis

Life Data Analysis

Introduction

Life Data Analysis (LDA), also known as Weibull Analysis or Reliability Analysis, is a statistical approach used to evaluate the time-to-failure data of products or systems. This method helps engineers and decision-makers understand how long a product is likely to function before failure, enabling informed predictions about reliability, maintainability, and performance. By fitting lifetime data to probability distributions---such as the Weibull, exponential, or lognormal---LDA provides critical metrics like mean time to failure (MTTF), failure rate, and reliability over time. These insights support better product design, warranty planning, and maintenance scheduling.

LDA plays a crucial role across industries where operational uptime and safety are essential, including aerospace, manufacturing, electronics, and energy. It allows organizations to make data-driven decisions by quantifying failure behavior and identifying the dominant failure modes. Additionally, LDA can incorporate both complete failure data and censored (incomplete) data, such as items still functioning at the end of a test or removed for other reasons. By doing so, it offers a flexible and powerful framework to improve system reliability and optimize lifecycle costs.

Statistical Background

Random Variable

A random variable is a numerical value that represents the outcome of a random experiment. It assigns numbers to the possible results of that experiment, allowing us to analyze and model uncertainty using mathematics. Random variables are fundamental in probability and statistics, as they help quantify outcomes and calculate probabilities.

There are two main types of random variables: discrete and continuous. A discrete random variable takes on countable values, such as 0, 1, 2, or 3 failures in a test, often arising from processes like counting or classification. In contrast, a continuous random variable can take on any value within a range or interval, such as the exact time to failure of a component. Continuous variables are typically measured rather than counted and are described using probability density functions (PDFs), while discrete variables use probability mass functions (PMFs).

In Life Data Analysis, we only handle continuous random variables.

Probability Density Function (PDF)

Assuming following is times-to-failure (days) data for a population of a 1000 Bearings.

Create a histogram with interval = 100 days as shown:

Suppose we want to estimate what percentage of the population would have failed by 300 days of operation—here's how we can do it:

There are 3 bars before 300 days. Its heights are 6, 31 and 74 respectively. Sum up the heights of these 3 bars, and we get 111 failures.

The total number of this sample is 1000 units. Therefore, the probability of a bearing failure is 0.111 for a 300 days operation.

Note that we can only query the probability of failures at 100, 200... at the resolution of the bin-size.

Let's overlay a function over the histogram as shown:

The function should have the following characteristics:

1. The function profile averages out the value of each bar.
2. The total area under this function is 1.

With this function, we can query the probability of failure at any time of interest. For example, you may query the probability of failure at 150 days by integrate the function from 0 to 150 days.

This function is called Probability Density Function. We denote pdf as f(t).

A Probability Density Function (PDF) is a mathematical function that describes the likelihood of a continuous random variable taking on a particular value. The probability that the variable falls within a specific range is found by calculating the area under the PDF curve over that interval. The total area under the entire PDF curve is always equal to 1, reflecting the certainty that the variable takes on some value within its domain.

Cumulative Density Function (CDF)

PDF provides a visualization of the distribution of the random variable, but trying read the area under the function up to a point in time can be a challenge.

CDF allows user to read the probability of failure at time t.

Most classical statistic books use F(t) to denote CDF. In reliability text book, Q(t) is used. CDF is also known as Unreliability function.

As shown in the figure below, we can read from CDF: Q(t=300 days) = 0.11

To read from the PDF, you must determine the area under the graph f(t) from t=0 to t=300 days. This is not a trivial task!

The CDF is the integral of the PDF.

$Q(t_{0}) = P(t \leq t_{0}) = \int_{- \infty}^{t_{0}}{f(t)dt}$

Conversely, The PDF is the derivative of the CDF:

$f(t) = \frac{d(Q(t))}{dt}$

Reliability, Unreliability

The unreliability plot (CDF) allows users to read the probability at time t directly.

To obtain the Reliability at time t,

$R(t) = 1 - Q(t)$

Using the equation, the reliability plot can be derived.

Failure Rate

Mathematically, failure rate at time t, λ(t), is defined as the probability of failure per unit time, given that the population have survived till time t.

From definition,

where N(t) is the number of survivals at t, and $N_{0}$ is the initial population size.

Condition Reliability

Conditional reliability $R(t/T)$is the probability of a unit surviving a mission of 𝑡 duration, given that it has already accumulated an age of 𝑇.

For an item to survival from point A to point C, it must have survived from A to B, and B to C. Hence,

Mean and Median

Mean time to failure (MTTF), $\overline{T}$ is defined as:

$\overline{T} = \int_{0}^{\infty}t \bullet f(t)dt$

Median is the time $\breve{T}$ by which 50% of the population fails.

$Q\left( \breve{T} \right) = 0.5$

Reliability Data Classification

In Life Data Analysis, reliability data falls into the following main types:

1. Complete (Failure) Data: Each item is observed until failure. The exact failure time is known. This is the most straightforward type of data.
2. Right-Censored Data (Suspension Data): The item did not fail during the observation period. We only know it survived up to a certain time. This is common in real-world testing or field data.
3. Interval-Censored Data: The exact failure time is unknown but falls within a known time interval (e.g., checked periodically and found failed between two inspections).
4. Left-Censored Data: The failure occurred before a known time, but the exact time is not observed.

These data types are used together to estimate life distributions and reliability metrics.

Distribution Models

Distribution models are used to statistically describe the time until failure of components, systems, or products. These models help engineers and analysts understand the underlying failure behavior---whether failures occur randomly, increase with age, or follow a particular degradation pattern. By fitting a suitable distribution to observed failure data, one can estimate key reliability metrics such as failure rate, reliability over time, and mean time to failure. This allows for informed decisions in design, maintenance, warranty forecasting, and risk management, ultimately improving product reliability and operational efficiency.

Exponential (λ)

The exponential distribution is commonly used to model the time between random, memoryless events---particularly when the failure rate is constant over time. While simple, it's important and often serves as a baseline model in reliability and survival analysis.

Probability Density Function, PDF

The exponential 𝑝𝑑𝑓 𝑓(𝑡) is given by:

$f(t) = {\lambda e}^{- \lambda t}$

If $\lambda$ is decreased, the Exponential distribution stretches to the right. The height decreases to preserve the total area under the probability density function (PDF), which must always equal 1.

Cumulative Density Function, CDF

If $\lambda$ is decreased, the Exponential distribution stretches to the right.

Failure Rate Function

The mean, $\overline{T}$, is given by:

Mean Time to Failure (MTTF)

The mean, $\overline{T}$, is given by:

Median Life

The median life, $\breve{T}$ is the time by which 50% of the population fails.

Weibull (β,η)

The Weibull distribution is commonly used to model time to failure data. Its flexibility allows it to represent a wide range of failure behaviors depending on its shape parameter β.

Probability Density Function, PDF

The Weibull PDF 𝑓(𝑡) is given by:

$f(t;\beta,\eta) = \frac{\beta}{\eta}\left( \frac{t}{\eta} \right)^{\beta - 1}e^{{- \left( \frac{t}{\eta} \right)}^{\beta}}$

Effect of $\beta$ and $\eta$ on the PDF:

For $0 < \beta < 1$: $f(t)$ is Infinity at start, decays rapidly

For $\beta = 1$: Exponential shape

For $\beta > 1$: $f(t) = 0$ at start, peaks then decline

If 𝜂 is increased while 𝛽 remains constant, the Weibull distribution stretches to the right. The overall shape of the curve stays the same, but the peak height decreases to preserve the total area under the probability density function (PDF), which must always equal 1.

Cumulative Density Function, CDF

The CDF, also known as unreliability, is commonly denoted as Q(t) or F(t).

$Q(t;\beta,\eta) = 1 - e^{- \left( \frac{t}{\eta} \right)^{\beta}}$

Effect of $\beta$ and $\eta$ on the CDF:

For $0 < \beta < 1$: Rapid rise early, levels off slowly

For $\beta = 1$: Exponential shape

For $\beta > 1$: Shaped curve (slow start, then rapid rise)

If 𝜂 is increased while 𝛽 remains constant, the CDF curve stretches to the right. The overall shape of the curve stays the same.

Failure Rate Function

The failure rate function $\lambda(t)$provides the percentage of failures occurring per unit time:

$\lambda(t) = \frac{f(t)\text{~}}{R(t)}$

Substituting $f(t) = \frac{\beta}{\eta}\left( \frac{t}{\eta} \right)^{\beta - 1}e^{{- \left( \frac{t}{\eta} \right)}^{\beta}}$

and $R(t) = e^{- \left( \frac{t}{\eta} \right)^{\beta}}$:

Effect of $\beta$ and $\eta$ on the failure rate function:

For $0 < \beta < 1$: High failure rate initially, and decrease thereafter. Ideal for modelling infant mortality

For $\beta = 1$: Constant failure rate

For $\beta > 1$: $f(t) = 0$ S-shaped curve (slow start, then rapid rise)

Decreasing $\eta$ shifts failures earlier, increasing the failure rate at the same time t.

Mean Time to Failure (MTTF)

The mean, $\overline{T}$, of the Weibull pdf is given by:

$\overline{T} = \eta \bullet \Gamma\left( \frac{1}{\beta} + 1 \right)$

Where

$\Gamma(n) = \int_{0}^{\infty}{e^{- x}x^{n - 1}dx}$

Median Life

The median life, $\breve{T}$ is the time by which 50% of the population fails.

Normal (μ,σ)

The normal distribution is a symmetric, bell-shaped probability distribution that is widely used in statistics. In life data analysis, it's used less frequently than Weibull or lognormal because it allows negative values, which are not realistic for time-to-failure. However, it can still be appropriate in certain cases---particularly when variability is low and failures are clustered around a central value.

Probability Density Function, PDF

The Normal PDF 𝑓(𝑡) is given by:

$f(t;\mu,\sigma) = \frac{1}{\sigma\sqrt{2\pi}}e^{- \frac{1}{2}\left( \frac{t - \mu}{\sigma} \right)^{2}}$

Effect of $\mu$ and $\sigma$ on the PDF:

$\mu$ determines the center (mean) or location of the distribution.

$\sigma$ controls the spread or dispersion of the data around$\ \mu$.

Cumulative Density Function, CDF

$Q(t) = \frac{1}{\sigma\sqrt{2\pi}}\int_{- \infty}^{t}e^{- \frac{1}{2}\left( \frac{x - \mu}{\sigma} \right)^{2}}\mathbb{d}x$

Alternatively, Normal CDF can also be expressed using the standard normal CDF:

$Q(t) = \Phi\left( \frac{t - \mu}{\sigma} \right)$

Where

$\Phi( \bullet )$: standard normal CDF

Effect of $\mu$ and $\sigma$ on the CDF:

Failure Rate Function

The lognormal failure rate is given by:

$\lambda(t) = \frac{f(t)}{R(t)}$

where

$f(t) = \frac{1}{\sigma\sqrt{2\pi}}e^{- \frac{1}{2}\left( \frac{t - \mu}{\sigma} \right)^{2}}$

$R(t) = 1 - \frac{1}{\sigma\sqrt{2\pi}}\int_{- \infty}^{t}e^{- \frac{1}{2}\left( \frac{x - \mu}{\sigma} \right)^{2}}\mathbb{d}x$

Mean Time to Failure (MTTF) and Median Life

Since the normal distribution is symmetrical, the median is equal to the mean, and is the distribution parameter $\mu$:

$\overline{T} = \breve{T} = \mu$

Lognormal (μ',σ')

The lognormal distribution is used to model failure times where the logarithm of the data follows a normal distribution. This distribution is particularly suitable for products or systems whose failure mechanisms are driven by the accumulation of damage or degradation over time---such as fatigue, corrosion, or chemical breakdown. The lognormal model produces a right-skewed curve, meaning most units fail around a certain time, but a few last significantly longer. It is characterized by a slow initial failure rate that increases over time, making it useful for modeling wear-out failure behavior where variability in time-to-failure is high.

Probability Density Function, PDF

The Lognormal PDF 𝑓(𝑡) is given by:

$f(t) = \frac{1}{{t \bullet \sigma}'\sqrt{2\pi}}e^{- \frac{1}{2}\left( \frac{ln(t) - \mu'}{\sigma'} \right)^{2}\text{~}}$

Where

$\mu'$: mean of the natural log of 𝑡 (log-location parameter)

$\sigma'$: standard deviation of the natural log of 𝑡 (log-scale parameter)

$\sigma'$ controls the spread and skewness of the distribution.

A larger $\sigma'$ value makes the curve flatter and wider, indicating a greater likelihood of very long lifetimes.

A smaller $\sigma'$ value (less than 1) results in a narrower and taller curve, meaning failures are more concentrated around the median.

The parameter $\mu'$ controls the horizontal positioning (time scaling) of the PDF but does not affect its shape.

Cumulative Density Function, CDF

$Q(t) = \frac{1}{t\sigma'\sqrt{2\pi}}\int_{0}^{t}{e^{- \frac{1}{2}\left( \frac{ln(s) - \mu'}{\sigma'} \right)^{2}\text{~}}ds}$

Where

$\mu'$: mean of the natural log of 𝑡 (log-location parameter)

$\sigma'$: standard deviation of the natural log of 𝑡 (log-scale parameter)

Alternatively, Lognormal CDF can also be expressed as:

$Q(t) = \Phi\left( \frac{\ln(t) - \mu'}{\sigma'} \right)$

Where

$\Phi( \bullet )$: standard normal CDF

A larger $\sigma'$ stretches the CDF horizontally, indicating that failures are spread out over a wider range of times.

A smaller $\sigma'$ compresses the CDF and makes it steeper, meaning that failures are more concentrated within a narrower time window.

The parameter $\mu'$ controls the horizontal positioning (time scaling) of the CDF but does not affect its shape.

Failure Rate Function

The lognormal failure rate is given by:

$\lambda(t) = \frac{f(t)}{R(t)}$

where

$f(t) = \frac{1}{{t \bullet \sigma}'\sqrt{2\pi}}e^{- \frac{1}{2}\left( \frac{ln(t) - \mu'}{\sigma'} \right)^{2}\text{~}}$

$R(t) = 1 - \frac{1}{t\sigma^{'\sqrt{2\pi}}}\int_{0}^{t}{e^{- \frac{1}{2}\left( \frac{\ln(s) - \mu'}{\sigma'} \right)^{2}\text{~}}ds}$

Mean Time to Failure (MTTF)

The mean, $\overline{T}$, is given by:

This integral evaluates to a known result in statistics:

$\overline{T} = e^{\mu + \frac{1}{2}{\sigma'}^{2}}$

Median Life

The median life, $\breve{T}$ is the time by which 50% of the population fails.

$Q\left( \breve{T} \right) = \Phi\left( \frac{\ln\left( \breve{T} \right) - \mu'}{\sigma'} \right) = 0.5$

Since $\Phi(0) = 0.5$,

Parameter Estimation

The goal of parameter estimation is to determine the values of distribution parameters that best characterize a given set of observed data.

The following sections describe two widely used methods for this purpose: Probability Plotting and Maximum Likelihood Estimation (MLE). These techniques are commonly applied to fit parametric models such as the Weibull, lognormal, normal, and exponential distributions using sample data.

Probability Plotting

Probability Plot is a graphical method to estimate the parameters (e.g., mean, standard deviation, shape) of a theoretical probability distribution that best fits a given dataset.

The following dataset is used for illustration. {30, 40, 20, 50}.

We first sort the Failure Data: {20, 30, 40, 50}.

For each failure time ti​, compute the cumulative probability qi using Cumulative Binomial equation.

$0.5 = \sum_{k = j}^{N}{\binom{N}{k}q^{k}(1 - q)^{N - k}}$

A cumulative probability $q_{j}$ is associated to each $t_{j}$.   The next step is to create the Probability paper. The X and Y axis are transformed such that the CDF appears as a straight line when plotted in this graph paper.

The following sections describe the concept of creating probability paper and solving for the various distribution (Exponential, Weibull, Normal and Lognormal) parameters.

Exponential Distribution

The cumulative distribution function (CDF) is:

Taking the natural log of both sides:

$\ln\left( R(t) \right) = - \lambda t$

Hence, the plot of $\ln\left( R(t) \right)$ against t with form a straight line with slop -λ.

The co-ordinates (t, ln(R(t)) are plotted on the Probability-Exponential paper as shown.

Since R(t=0) is 1, the straight is forced to pass through

$(t = 0,\ ln\left( R(0) \right) = 0)$ as shown.

To obtain the slop, pick a point on the fitted line know Y value to calculate the slop.
At Y = -2.303, the corresponding time is t = 75.

The Exponential parameter is:

$\lambda = - slop = - \frac{0 - \left( - \text{2.303~} \right)}{75} = 0.0307$ per unit time

Weibull Distribution

For the Weibull distribution, the cumulative density function can be written as:

$\Rightarrow \ ln\left( \ln\left( \frac{1}{1 - Q(t)} \right) \right) = \beta \bullet ln(t) - \beta \bullet ln(\eta)$

Linearize the equation: $Y = \beta X + C$,

where

$Y = ln\left( \ln\left( \frac{1}{1 - Q(t)} \right) \right)$

$X = ln(t)$

$slop = \beta$

Hence, the plot of Y against X with form a straight line with slop $\beta$.

The co-ordinates (X, Y) are plotted on the Probability-Weibull paper as shown.

From Weibull CDF,

$Q(t) = 1 - e^{- \left( \frac{t}{\eta} \right)^{\beta}}$

Let $t = \eta$

$Q(\eta) = 1 - e^{- (1)^{\beta}} = 63.1\%$

Therefore, the time corresponding to Q=63.1% is the η value (=40), as shown in the graph below.

To obtain σ, pick two points on the fitted line with known Y-values (e.g., Y=0.834 and Y=-0.367), and calculate the slop.

For Z=0.834, cumulative probability Q(t)=0.9

For Z=-0.367, cumulative probability Q(t)=0.5

The $\beta$ is:

$\beta = \frac{0.834 - ( - 0.367)}{\ln(55.4) - ln(34.6)} = 2.56$

Normal Distribution

For the normal distribution $Nor(\mu,\sigma)$, the cumulative density function can be written as:

where $\Phi^{- 1}$ is the inverse of the standard normal CDF.

If we set $\Phi^{- 1}\left( Q(t) \right)$ as the Y-axis, and t (no transformation) as the X-axis, then the CDF of a Normal distribution will become a straight line with y-intercept = $- \frac{\mu}{\sigma}$ and slop =$\frac{1}{\sigma}$ on this transform graph paper.

Now, we transform the cumulative probabilities Q(t) into standard normal quantiles (Z-values):
$Z = \Phi^{- 1}\left( Q(t) \right)$

Note: You can use a standard normal table or a function like NORM.S.INV(Q) in Excel to get the Z values.

The co-ordinates (t, Z) are plotted on the Probability-Normal paper as shown.

Mathematically, we can derive the mean and standard-deviation from the y-interception ($= - \frac{\mu}{\sigma}$), and the slop ($\frac{1}{\sigma}$).

For Normal distribution, $Q(t = \mu) = 50\%$. From the graph, $\mu = 35$.

To obtain σ, pick two points on the fitted line with known Z-values (e.g., Z=1 and Z=-1 ) and calculate the inverse of the slop.

For Z=1, cumulative probability Q(t)=0.841

For Z=-1, cumulative probability Q(t)=0.159

Note: You can use a standard normal table or a function like NORM.S.DIST(Z) in Excel to get the cumulative probability Q(t) values.

The Standard-deviation is:

Lognormal Distribution

If a variable t follows a lognormal distribution, then ln(t) follows a Normal distribution.

So, we simply take logarithms of the failure times and then plot them against standard normal quantiles.

For the Lognormal distribution $Log(\mu',\sigma')$, the cumulative density function can be written as:

$Q(t) = \Phi\left( \frac{\ln(t) - \mu'}{\sigma'} \right)$

$\Rightarrow \ \Phi^{- 1}\left\lbrack Q(t) \right\rbrack = - \frac{\mu'}{\sigma'} + \frac{1}{\sigma'}\ln(t)$

where $\Phi^{- 1}$ is the inverse of the standard normal CDF.

$\mu'$ is the Mean of the log-transformed failure times.

$\sigma'$ is the Standard deviation of the log-transformed failure times.

If we set $\Phi^{- 1}\left( Q(t) \right)$ as the Y-axis, and ln(t) as the X-axis, then the CDF of a Lognormal distribution will become a straight line with a slope $= \frac{1}{\sigma'}$

Now, we transform the cumulative probabilities Q(t) on the Y-axis into standard normal quantiles (Z-values): $Z = \Phi^{- 1}\left( Q(t) \right)$, and the X-axis into ln(t).

Note: You can use a standard normal table or a function like NORM.S.INV(Q) in Excel to evaluate $\Phi^{- 1}\left( Q(t) \right)$ the Z value in the table below.

The co-ordinates (ln(t), Z) are plotted on the Probability-Lognormal paper as shown.

The log-mean (µ') occurs at Z=0, i.e., Q(t)=50%. At this point $\ln(t) = \mu$

The corresponding t for Q(t)=50% is 33.1, hence $\mu\ = \ln(33.1) = 3.5$

To obtain σ', pick two points on the fitted line with known Z-values (e.g., Z=1 and Z=-1), and calculate the reciprocal of the slop.

For Z=1, cumulative probability Q(t)=0.841

For Z=-1, cumulative probability Q(t)=0.159

Note: You can use a standard normal table or a function like NORM.S.DIST(Z) in Excel to get the cumulative probability Q(t) values.

The corresponding t values at those points are 18.3 and 52.5 as shown.

The Log Standard-deviation is:

Rank Regression

by fitting it to the observed data. The goal is to find the line that minimizes the sum of squared differences between the observed values and the values predicted by the line.

Rank Regression on Y

The RRY line (Rank Regression on Y) is a regression line that minimizes the horizontal distance (errors in the Y-direction) between each data point and the regression line.

It is used when we want to predict the time (e.g., failure time) for a given probability (e.g., 90% unreliability).

Given paired data points: (x1, y1), (x2, y2) ..., (xn, yn), and that the x-values are known exactly.

$y_{i} = a' + b'x_{i} + e_{i}$

Where:

$b' = \frac{\sum_{i = 1}^{n}{\left( y_{i} - \overline{y} \right)\left( x_{i} - \overline{x} \right)}}{\sum_{i = 1}^{n}\left( x_{i} - \overline{x} \right)^{2}}$

$a' = \overline{y} - b^{'\overline{x}}$

$\overline{x} = \frac{1}{n}\sum_{i = 1}^{n}x_{i}$

$\overline{y} = \frac{1}{n}\sum_{i = 1}^{n}y_{i}$

$e_{i}$ is the error term.

Rank Regression on X

The RRX line (Rank Regression on X) is a regression line that minimizes the horizontal distance (errors in the X-direction) between each data point and the regression line.

It is used when we want to predict the probability (e.g., 90% unreliability) for a given time.

Given paired data points: (x1, y1), (x2, y2) ..., (xn, yn), and that the y-values are known exactly.

To minimize the error in x-direction, regression of x on y is performed, i.e.,

$x_{i} = a' + b'y_{i} + e_{i}$

Where:

$b' = \frac{\sum_{i = 1}^{n}{\left( y_{i} - \overline{y} \right)\left( x_{i} - \overline{x} \right)}}{\sum_{i = 1}^{n}\left( y_{i} - \overline{y} \right)^{2}}$

$a' = \overline{x} - b'\overline{y}$

$\overline{x} = \frac{1}{n}\sum_{i = 1}^{n}x_{i}$

$\overline{y} = \frac{1}{n}\sum_{i = 1}^{n}y_{i}$

$e_{i}$ is the error term.

Correlation Coefficient

One key concept in rank regression is the correlation coefficient, which serves as a measure of how well the chosen distribution fits the data.

The correlation coefficient ρ is defined as:

Where:

• $\overline{x}$ and $\overline{y}$​ are the sample means of x and y,
• $Cov(X,Y)$ is the covariance,
• $\sigma_{X}$​, $\sigma_{Y}$ are the standard deviations of x and y.

Interpretation

• The range is between -1 and +1.
• -1 is a perfect fit with negative slop.
• +1 is a perfect fit with positive slop.

Consider the following dataset:

{20, 30, 40, 50}.

The following plot shows the corresponding $\rho$ values when fitted to Exponential, Weibull, Normal and Lognormal model using Rank Regression on X (RRX).

In the case of probability plot method (rank regression), the $\rho$ value can be used as an indicator of how well the chosen distribution fits the data.

Maximum Likelihood Estimations

Maximum Likelihood Estimation (MLE) is a statistical method used to estimate the parameters of a probability distribution by finding the values that make the observed data most probable under the assumed distribution. The key idea is to construct a likelihood function, which represents the probability (or probability density) of the observed data as a function of the unknown parameters. MLE seeks the parameter values that maximize this likelihood function.

Let:

$f(t;\theta)$: Probability density function (PDF) of the distribution (e.g., Weibull, lognormal).

$F(t;\theta)$: Cumulative distribution function (CDF).

$R(t;\theta) = 1 - F(t;\theta)$: Reliability function (survival probability).

$\theta$: Vector of unknown distribution parameters (e.g., shape and scale).

$n_{1}$: exact failures at times $t_{i}$.

$n_{2}$: right-censored observations at times $t_{i}$.

$n_{3}$: interval-censored observations between $t_{Li}$ and $t_{Lu}$.

The complete likelihood function is given by:

$L(\theta) = \prod_{i = 1}^{n_{1}}{f\left( t_{i};\theta \right) \times \text{~}}\prod_{i = 1}^{n_{2}}{R\left( t_{i};\theta \right) \times \prod_{i = 1}^{n_{3}}{\left\lbrack F\left( t_{U_{i}};\theta \right) - F\left( t_{L_{i}};\theta \right) \right\rbrack\text{~}}}$

Since maximizing the product is mathematically complex, it's easier to work with the log-likelihood:

$\Lambda = ln\left( L(\theta) \right) = \sum_{i = 1}^{n_{1}}{\ln\left( f\left( t_{i};\theta \right) \right) \times \sum_{i = 1}^{n_{2}}{\ln\left( R\left( t_{i};\theta \right) \right) \times \sum_{i = 1}^{n_{3}}{\ln\left\lbrack F\left( t_{U_{i}};\theta \right) - F\left( t_{L_{i}};\theta \right) \right\rbrack\text{~}}}}$

MLE finds the parameter values θ that maximize this log-likelihood,Λ.

Let's consider fitting the following dataset to Weibull distribution using MLE analysis method:

PDF:

$f(t;\beta,\eta) = \frac{\beta}{\eta}\left( \frac{t}{\eta} \right)^{\beta - 1}e^{{- \left( \frac{t}{\eta} \right)}^{\beta}}$

CDF:

$Q(t;\beta,\eta) = 1 - e^{- \left( \frac{t}{\eta} \right)^{\beta}}$

Reliability:

$R(t;\beta,\eta) = e^{- \left( \frac{t}{\eta} \right)^{\beta}}$

Log-Likelihood Function:

$\small\Lambda = ln\left( f(20;\beta,\eta) \right) + ln\left( f(30;\beta,\eta) \right) + 2 \times ln\left( R(40;\beta,\eta) \right) + ln\left\lbrack F(30;\beta,\eta) - F(50;\beta,\eta) \right\rbrack$

Substituting Weibull functions:

$\small\Lambda = \left\lbrack \ln\left( \frac{\beta}{\eta} \right) + (\beta - 1)\ln(20) - \left( \frac{20}{\eta} \right)^{\beta} \right\rbrack + \left\lbrack \ln\left( \frac{\beta}{\eta} \right) + (\beta - 1)\ln(30) - \left( \frac{30}{\eta} \right)^{\beta} \right\rbrack$ $- 2\left( \frac{40}{\eta} \right)^{\beta} + ln\left( e^{- \left( \frac{30}{\eta} \right)^{\beta}} - e^{- \left( \frac{50}{\eta} \right)^{\beta}} \right)$

The above Log-Likelihood is a function of $\beta,\eta$.

Maximum Likelihood Estimation (MLE) determines the parameter values $\beta$ and $\eta$ that maximize the log-likelihood function $\Lambda$.

Because this function is often nonlinear and complex, numerical optimization methods are required. The Weibull-Toolbox uses the Nelder--Mead optimization method to perform this estimation.

For this example, the β and η values are 3.37 and 41.9 respectively.

If these values are substituted back into the Log-Likelihood equation, we get:

$\small\Lambda = \left\lbrack \ln\left( \frac{\text{3.37}}{\text{41.9}} \right) + \left( \text{3.37} - 1 \right)\ln(20) - \left( \frac{20}{\text{41.9}} \right)^{\text{3.37}} \right\rbrack + \left\lbrack \ln\left( \frac{\text{3.37}}{\text{41.9}} \right) + \left( \text{3.37} - 1 \right)\ln(30) - \left( \frac{30}{\text{41.9}} \right)^{\text{3.37}} \right\rbrack$ $- 2\left( \frac{40}{\text{41.9}} \right)^{\text{3.37}}+ ln\left( e^{- \left( \frac{30}{\text{41.9}} \right)^{\text{3.37}}} - e^{- \left( \frac{50}{\text{41.9}} \right)^{\text{3.37}}} \right) {= - 10.28}$

In Weibull-Toolbox, the log-likelihood value is displayed under Analysis Results as the LK Value.

Suppose we fit the same dataset to several distributions---Normal, Lognormal, Weibull, and Exponential---and compare their LK Values:

Since the Normal distribution yields the highest LK Value (i.e., the least negative), it represents the best fit among the options. In contrast, the Exponential distribution shows a significantly lower LK Value, indicating a poor fit to the data.

Confidence Bounds

Selecting a method. No single confidence-bound method is universally preferred. Fisher Matrix, likelihood-ratio, and simulation methods use different approximations and may produce different results, particularly when data are limited or heavily censored. The selected method should be appropriate for the data, fitted model, and engineering question.

Fisher Matrix Confidence Bounds

After fitting a Weibull model to failure dataset, we can quantify the uncertainty in the parameter estimates ($\widehat{\beta}$, $\widehat{\eta}$) using Fisher Matrix Confidence Bounds (FMCBs) from which the variance for β and η can be derived.

Fisher Matrix Confidence Bounds on Parameter

We will demonstrate Fisher Matrix Confidence Bounds with Weibull distribution. The same concept is applicable to the other life distribution models.

The likelihood function for the 2P-Weibull distribution is:

$L(\beta,\eta) = \prod_{i = 1}^{n}\left\lbrack \left( \frac{\beta}{\eta}\left( \frac{t_{i}}{\eta} \right)^{\beta - 1} \right)^{\delta_{i}} \cdot exp\left( - \left( \frac{t_{i}}{\eta} \right)^{\beta} \right) \right\rbrack$

• $t_{i}$​ is the observed time for the i-th item,
• $\delta_{i}$ =1 for failures, and $\delta_{i}$ = 0 for censored data.

The log-likelihood function is:

$\ln\left( L(\beta,\eta) \right) = \sum_{i = 1}^{n}\delta_{i}\left\lbrack \ln(\beta) - \beta ln(\eta) + (\beta - 1)\ln\left( t_{i} \right) \right\rbrack - \sum_{i = 1}^{n}\left( \frac{t_{i}}{\eta} \right)^{\beta}$

Fisher Information Matrix for Weibull distribution

The Fisher Information Matrix I(β,η) quantifies the curvature of the log-likelihood function around the MLE point. It is defined as the negative expectation of the second derivatives:

Its inverse gives the covariance matrix of the estimates:

We assume that $ln(\beta)$ and $ln(\eta)$ are approximately normally distributed. It ensures that bounds for β and η remain strictly positive.

Shape parameter $\beta$ bounds:

$\ln\left( \beta_{U,L} \right) = ln\left( \widehat{\beta} \right) \pm z_{\frac{\alpha}{2}} \cdot \sqrt{Var\left( \ln\left( \widehat{\beta} \right) \right)}$

$Var(ln(\widehat{\beta})$ can be approximated using the Delta Method as:

$Var(ln(\widehat{\beta}) \approx \frac{Var(\widehat{\beta})}{\beta^{2}}$

Hence, $\ln\left( \beta_{U,L} \right) = \ln\left( \widehat{\beta} \right) \pm z_{\frac{\alpha}{2}} \cdot \frac{\sqrt{Var\left( \widehat{\beta} \right)}}{\widehat{\beta}}$

For 90% confidence level, α=0.1, and $z_{\frac{\alpha}{2}} = 1.64$

The $\beta$ 2-sided bound at confidence level α:

$\beta_{U,L} = \widehat{\beta} \bullet e^{\pm z_{\frac{\alpha}{2}}\frac{\sqrt{Var\left( \widehat{\beta} \right)}}{\widehat{\beta}}}$

The scale parameter $\eta$ bounds:

$\ln\left( \eta_{U,L} \right) = ln\left( \widehat{\eta} \right) \pm z_{\frac{\alpha}{2}} \cdot \frac{\sqrt{Var\left( \widehat{\eta} \right)}}{\widehat{\eta}}$

The $\eta$ 2-sided bound at confidence level $\alpha$:

$\eta_{U,L} = \widehat{\eta} \bullet e^{\pm z_{\frac{\alpha}{2}}\frac{\sqrt{Var\left( \widehat{\eta} \right)}}{\widehat{\eta}}}$

Example:

A dataset {20, 30 ,40, 50} is fitted to a 2-P Weibull distribution using RRX. Find the 90% confidence bounds on the parameters.

The RRX solution: {beta=2.55, eta=39.9}

Fisher Information Matrix:

Covariance Matrix:

$Var\left( \widehat{\beta} \right) = 1.215$

$Var\left( \widehat{\eta} \right) = 76.42$

$\beta_{U} = \widehat{\beta} \bullet e^{z_{\frac{\alpha}{2}}\frac{\sqrt{Var\left( \widehat{\beta} \right)}}{\widehat{\beta}}} = 2.55 \times e^{1.64\left( \frac{\sqrt{1.215}}{2.55} \right)} = 1.253$27.8

$\beta_{L} = \widehat{\beta} \bullet e^{- z_{\frac{\alpha}{2}}\frac{\sqrt{Var\left( \widehat{\beta} \right)}}{\widehat{\beta}}} = 2.55 \times e^{- 1.64\left( \frac{\sqrt{1.215}}{2.55} \right)} = 5.188$

$\eta_{U} = \widehat{\eta} \bullet e^{z_{\frac{\alpha}{2}}\frac{\sqrt{Var\left( \widehat{\eta} \right)}}{\widehat{\eta}}} = 39.9 \times e^{1.64\left( \frac{\sqrt{76.42}}{39.9} \right)} = 57.23$

$\eta_{L} = \widehat{\eta} \bullet e^{- z_{\frac{\alpha}{2}}\frac{\sqrt{Var\left( \widehat{\eta} \right)}}{\widehat{\eta}}} = 39.9 \times e^{- 1.64\left( \frac{\sqrt{76.42}}{39.9} \right)} = 27.84$

The concept of FM bound on an estimate can be summarized as follow:

1- Obtain the estimate through statistical analysis (In our case, MLE or RRX) and serve as our mean $\eta$.

2- Derive the standard deviation σ, from Fisher Information Matrix.

3- Bound the estimate by the normal distribution $N(\eta, \sigma)$ with a standard normal quantile defining the confidence level.

Fisher Matrix Confidence Bounds on Time

Consider a dataset {20, 30 ,40, 50} that is fitted to a 2-P Weibull distribution. Find the 90% confidence bounds for the time interval where 10% of the population fails (unreliability Q=10%).

The RRX solution: {$\beta = 2.55$, $\eta = 39.9$}

Weibull CDF:

$Q(t) = 1 - e^{- \left( \frac{t}{\eta} \right)^{\beta}}$

For $Q = 10\%$, $t = 16.5$

We bound a Normal distribution on time with mean = 16.5 as shown.

Since the x-axis is in natural-log scale Ln(t), the mean is actually ln(16.5).

We also need the variance Var(ln($\widehat{t}$)) to completely define the Normal bound.

Let u =ln(t), the upper and lower bounds on u:

$u_{U} = \widehat{u} + z_{\frac{\alpha}{2}} \cdot \sqrt{Var\left( \widehat{u} \right)}$

$u_{L} = \widehat{u} - z_{\frac{\alpha}{2}} \cdot \sqrt{Var\left( \widehat{u} \right)}$

We have Var($\widehat{\beta}$) and Var($\widehat{\eta}$), but not Var($\widehat{u}$). To derive Var($\widehat{u}$), we use Delta Method (Ref: Statistical Methods for Reliability Data by Meeker & Escobar, Appendix B, B2-Statistical Error Propagation-The Delta Method).

First, we need a relationship that expresses u in-terms of β and η.

From Weibull CDF

$Q(t) = 1 - e^{- \left( \frac{t}{\eta} \right)^{\beta}}$

Linearize the equation:

$lnln\left( \frac{1}{1 - Q(t)} \right) = \beta\left\lbrack ln(t) - ln(\eta) \right\rbrack$

$ln(t) = \frac{1}{\beta}lnln\left( \frac{1}{1 - Q(t)} \right) + ln(\eta)$

Substitute $u=ln(t)$

$u = \frac{1}{\beta}lnln\left( \frac{1}{1 - Q(t)} \right) + ln(\eta)$

Using Delta Method:

$Var(u) = \left( \frac{\partial u}{\partial\beta} \right)^{2}Var(\beta) + \left( \frac{\partial u}{\partial\eta} \right)^{2}Var(\eta) + 2\left( \frac{\partial u}{\partial\beta} \right)$$\left( \frac{\partial u}{\partial\eta} \right)$ $\small Cov(\beta,\eta$)

$= \frac{1}{\beta^{4}}\left( lnln\left( \frac{1}{(1 - Q)} \right) \right)^{2}Var(\beta) + \left( \frac{1}{\eta} \right)^{2}Var(\eta) - \frac{2}{{\eta\beta}^{2}}\left( lnln\left( \frac{1}{(1 - Q)} \right) \right)Cov(\beta,\eta)$

$Var\left( \widehat{\beta} \right) = 1.215$

$Var\left( \widehat{\eta} \right) = 76.42$

$Q = 0.1$

$\widehat{\beta} = 2.55$

$\widehat{\eta} = 39.9$

$\therefore Var\left( \widehat{u} \right) = 0.1854$

$\widehat{u} = {ln(t)}{=ln(16.5)}{= 2.804}$

Now, we can apply:

$u_{U} = \widehat{u} + z_{\frac{\alpha}{2}} \cdot \sqrt{Var\left( \widehat{u} \right)} = 3.513$

$u_{L} = \widehat{u} - z_{\frac{\alpha}{2}} \cdot \sqrt{Var\left( \widehat{u} \right)} = 2.096$

Since $u= ln(t)$ => $t=exp(u)$

$t_{U}$ = exp($u_{U}$) = 33.5

$t_{L}$= exp($u_{L}$) = 8.13

The 90% confidence interval on t is evaluated at Q=10% We can repeat for different Q value:

By joining the points on the probability plot, we obtain the confidence bounds as shown.

Fisher Matrix Confidence Bounds on Reliability

A dataset {20, 30 ,40, 50} is fitted to a 2-P Weibull distribution. Find the 90% confidence bounds for unreliability at time t = 20.

The RRX solution: {beta=2.55, eta=39.9}

From Weibull CDF:

$Q(t) = 1 - e^{- \left( \frac{t}{\eta} \right)^{\beta}}$

For t=20, Q= 15.8%

We bound a Normal distribution on unreliability (y-axis) with mean = 15.8% as shown.

Since the y-axis is in double-log-reciprocal scale ($lnln\left( \frac{1}{1 - Q(t)} \right)$), the mean is actually $lnln\left( \frac{1}{1 - 0.158} \right)$.

We also need the variance $Var\left( lnln\left( \frac{1}{1 - \widehat{Q}(t)} \right) \right)$ to completely define the Normal bound.

Let

$u = lnln\left( \frac{1}{1 - Q(t)} \right)$, the upper and lower bounds on u:

$u_{U} = \widehat{u} + z_{\frac{\alpha}{2}} \cdot \sqrt{Var\left( \widehat{u} \right)}$

$u_{L} = \widehat{u} - z_{\frac{\alpha}{2}} \cdot \sqrt{Var\left( \widehat{u} \right)}$

We need a relationship that expresses $u$ in-terms of β and 𝜂, so that we can apply Delta Method to derive $Var(u)$.

From Weibull CDF:

$Q(t) = 1 - e^{- \left( \frac{t}{\eta} \right)^{\beta}}$

Linearize the equation:

$lnln\left( \frac{1}{1 - Q(t)} \right) = \beta\left\lbrack ln(t) - ln(\eta) \right\rbrack$

Substitute $\ u = lnln\left( \frac{1}{1 - Q(t)} \right)$

$u = \beta\left\lbrack ln(t) - ln(\eta) \right\rbrack$

Applying Delta Method:

$Var(u) = \left( \frac{\partial u}{\partial\beta} \right)^{2}Var(\beta) + \left( \frac{\partial u}{\partial\eta} \right)^{2}Var(\eta) + 2\left( \frac{\partial u}{\partial\beta} \right)$$\left( \frac{\partial u}{\partial\eta} \right)$ $\small Cov(\beta,\eta)$

$= \left( \ln(t) - ln(\eta) \right)^{2}Var(\beta) + \left( \frac{\beta}{\eta} \right)^{2}Var(\eta) - \frac{2\beta}{\eta}\left( \ln(t) - ln(\eta) \right)Cov(\beta,\eta)$

$Var\left( \widehat{\beta} \right) = 1.215$

$Var\left( \widehat{\eta} \right) = 76.42$

$t = 20$

$\widehat{\beta} = 2.55$

$\widehat{\eta} = 39.9$

$\therefore Var\left( \widehat{u} \right) = 0.8513$

$\widehat{u} = \beta{\left\lbrack ln(t) - ln(\eta) \right\rbrack = - 1.762}$

Now, we can apply:

$\small u_{U} = \widehat{u} + z_{\frac{\alpha}{2}} \cdot \sqrt{Var\left( \widehat{u} \right)} = - 0.2447$

$\small u_{L} = \widehat{u} - z_{\frac{\alpha}{2}} \cdot \sqrt{Var\left( \widehat{u} \right)} = - 3.280$

Since $u = lnln\left( \frac{1}{1 - Q(t)} \right)$ => $Q(t) = \ e^{{- e}^{u}}$

$\small q_{U} = \ e^{{- e}^{u_{U}}} = 0.543$

$\small q_{L} = \ e^{{- e}^{u_{L}}} = 0.0369$

The 90% confidence interval on Q is evaluated at t=20. We can repeat for different t values, and overlay it on the probability plot.

Likelihood Ratio Confidence Bounds

The Likelihood Ratio Principle

The Likelihood Ratio Confidence Bounds method is based on the likelihood ratio test.

The objective to estimate the confidence bounds on the β of Weibull distribution (that your dataset was fitted to).

Mathematically, the likelihood ratio statistic is:

$\small\Lambda(\beta) = - 2\ln\left( \frac{L\left( \beta,\widehat{\eta} \right)}{L\left( \widehat{\beta},\widehat{\eta} \right)} \right) = - 2\left\lbrack \ln\left( L\left( \beta,\widehat{\eta} \right) \right) - ln\left( L\left( \widehat{\beta},\widehat{\eta} \right) \right) \right\rbrack$

$L\left( \beta,\widehat{\eta} \right)$: the likelihood function with eta value fixed at $\widehat{\eta}$.

$L\left( \widehat{\beta},\widehat{\eta} \right)$: the likelihood value evaluated at $\widehat{\beta},\widehat{\eta}$.

The likelihood function for the 2P-Weibull distribution is:

$\small L(\beta,\eta) = \prod_{i = 1}^{n}\left\lbrack \left( \frac{\beta}{\eta}\left( \frac{t_{i}}{\eta} \right)^{\beta - 1} \right)^{\delta_{i}} \cdot exp\left( - \left( \frac{t_{i}}{\eta} \right)^{\beta} \right) \right\rbrack$