Random Variable Generation from Failure Rates

Md Aktar Ul Karim
Statistical Simulation Course
Symbiosis Statistical Institute, Pune

Contents

1 Introduction 1

2 Exponential Distribution for Constant Failure Rate: 2

3 Weibull Distribution for Variable Failure Rates: 3

4 Log-Normal Distribution: Non-Constant Failure Rate: 5

5 Exercises: 7

1. Introduction

Now, we will generate random variables from failure rates, which is essential in reliability analysis, survival
studies, and various applications in engineering and biological systems. In this case, we will focus on three key
distributions: Exponential, Weibull, and Log-Normal distributions. First, we will discuss their properties, and
then how can we model the failure rates, and then how to generate random variables using these distributions
in R.
What is a Failure Rate?

• The failure rate describes the instantaneous rate of failure for a system/component at any given time,
often denoted as λ(t).

• Commonly used in reliability engineering, survival analysis, and risk assessment.

• Failure rate can be constant or variable depending on the system’s characteristics.

Types of Failure Rate:

• Constant Failure Rate: The probability of failure does not change over time. Generally modeled by
the Exponential distribution.
 • Increasing/Decreasing Failure Rate: Generally modeled by the Weibull distribution, where the
failure rate can change over time. Failure becomes more likely as time progresses for increasing cases and
failure becomes less likely over time for decreasing cases.

• Non-Constant Multiplicative Failure Rate: In this case, failure processes are driven by multiple
influencing factors and generally modeled by the Log-Normal distribution.

2. Exponential Distribution for Constant Failure Rate:

The Exponential distribution is commonly used to model the time between failures in systems with a
constant failure rate.
Probability Density Function (PDF) is given by

f (t) = λe−λt ,

where λ is the failure rate (also known as the rate parameter) and t is the time until failure.
Cumulative Distribution Function (CDF): F (t) = 1 − e−λt .
Properties

• The probability of failure in the next instant is independent of how long the system has been running.

• Mean and Variance: The mean time to failure is 1

λ and the variance is 1
λ2
.

Generating Exponentially Distributed Random Variables:

We can generate an exponentially distributed random variable by using the R/python packages available. We
may also use the Inverse Transform Sampling to generate:

1
T = − ln(U ),
λ

where U is a uniform random variable on [0, 1].

Examples:
Lightbulb lifetimes follow an exponential distribution with a failure rate of λ = 0.02 failures per hour.

• Simulate 1000 random lightbulb lifetimes and plot the histograms.

• Plot the first 100 lifetimes.

• Add a red dashed line at the theoretical mean lifetime.

Ans:

2
# Set the failure rate
lambda = 0.02

# Generate 1000 random lifetimes of lightbulbs

[Link](123)
lifetimes = rexp(1000, rate = lambda)

# Histogram plot of the lifetimes

hist(lifetimes, breaks = 30, probability = TRUE, col = "grey",
main = "Histogram of Lightbulb Lifetimes with Exponential PDF",
xlab = "Lifetime (hours)", ylab = "Density", border = "white")

# Add the theoretical exponential distribution curve

curve(dexp(x, rate = lambda), col = "red", lwd = 2, add = TRUE)

# Generate a simple random plot of 100 lifetimes

plot(1:100, lifetimes[1:100], type = "b", pch = 19, col = "blue",
main = "Random Lifetimes of 100 Lightbulbs",
xlab = "Lightbulb Index", ylab = "Lifetime (hours)")

# Add the theoretical mean lifetime line

abline(h = 1/lambda, col = "red", lwd = 2, lty = 2) # Mean lifetime line

Figure 1: Lightbulb Problem

3. Weibull Distribution for Variable Failure Rates:

The Weibull distribution is more flexible than the Exponential distribution, capable of modeling systems
with increasing or decreasing failure rates.
Probability Distribution Function (PDF):
 β−1  β
β t − t
f (t) = e η , t ≥ 0.
η η
Depending the values of the shape parameter β, failure rate may increase or decrease:

3
 • β < 1: Decreasing failure rate.

• β = 1: Constant failure rate (Exponential distribution).

• β > 1: Increasing failure rate.

Scale Parameter η determines the spread of the distribution.

The cumulative distribution function (CDF) of the Weibull distribution is given by:
 β
t
−
F (t) = 1 − e η
, t≥0

• For t < 0, F (t) = 0.

• As t → ∞, F (t) → 1.

Generating Weibull Random Variables:

We can generate Weibull-distributed random variables directly by the R packages and also by using the inverse
transformation method:
T = η(− ln(U ))1/β ,

where U is a uniform random variable on [0, 1].

Example:
Generate 1000 random variables for a system with an increasing failure rate where β = 2 and η = 5. Visualize
the results.

Figure 2: Weibull Distribution

# Set parameters
beta = 2 # Shape parameter
eta = 5 # Scale parameter

4
# Generate 1000 random variables
[Link](123) # For reproducibility
random_vars = rweibull(1000, shape = beta, scale = eta)

# Histogram plot
hist(random_vars, breaks = 20, probability = TRUE, col = "grey",
main = "Weibull Distributed Random Variables",
xlab = "Time (t)", ylab = "Density", border = "white")

# Theoretical Weibull PDF curve

curve(dweibull(x, shape = beta, scale = eta), col = "red", lwd = 2, add = TRUE)

# Generate 1000 uniform random variables between 0 and 1

[Link](123)
U = runif(1000, min = 0, max = 1)

# Use the inverse transform method to generate Weibull-distributed random variables

weibull_vars = eta * (-log(1 - U))^(1 / beta)

# histogram Plot
hist(weibull_vars, breaks = 20, probability = TRUE, col = "blue",
main = "Weibull Random Variables (Using Uniform Distribution)",
xlab = "Time (t)", ylab = "Density", border = "white")

# Theoretical Weibull PDF curve

curve(dweibull(x, shape = beta, scale = eta), col = "red", lwd = 2, add = TRUE)

4. Log-Normal Distribution: Non-Constant Failure Rate:

The Log-Normal distribution is used when the failure process is influenced by multiplicative factors, common
in biological and financial models.
Probability Density Function (PDF):

1 (ln t−µ)2
f (t) = √ e− 2σ 2 ,
tσ 2π

µ and σ are the mean and standard deviation of the underlying normal distribution. This distribution is useful
in situations where the failure rate increases exponentially or where the natural logarithm of the time to failure
is normally distributed.
The mean and variance of the log-normal distribution os given by:
σ2
• Mean: E[T ] = eµ+ 2

2 2
• Variance: Var(T ) = (eσ − 1)e2µ+σ

Properties:

5
 • Skewed Distribution: The Log-Normal distribution is right-skewed, meaning most values are clustered
near the lower bound.

• Mean and Variance: The mean and variance will determine the shape and spread of the distribution.

The Cumulative Distribution Function (CDF) of the log-normal distribution is given by:
 
ln(t) − µ
F (t) = Φ , t > 0,
σ

Where Φ(·) is the CDF of the standard normal distribution.

Generating Log-Normal Random Variables: For this also we can generate random variables by using
packages or by using the inverse transformation algorithms:

−1 (U )
T = eµ+σΦ ,

where U is a uniform random variable on [0, 1].

Applications:

• Modeling the lifespan of animals in a population where survival depends on multiplicative environmental
and genetic factors.

• For some electronic components, we may experience failure rates that vary due to multiple factors, such
as environmental conditions and manufacturing differences.

Log-normal distribution models such cases.

Example Simulate 1000 random lifetimes using a Log-Normal distribution where µ = 1.5 and σ = 0.3.
Plot and analyze the distribution.

Figure 3: Log-Normal Distribution

6
# Parameters set up
mu = 1.5 # Mean of the log
sigma = 0.3 # Standard deviation of the log

# Simulate 1000 random lifetimes

[Link](123) # For reproducibility
random_lifetimes = rlnorm(1000, meanlog = mu, sdlog = sigma)

# Histogram plot
hist(random_lifetimes, breaks = 20, probability = TRUE,
col = "grey", main = "Histogram of Simulated Log-Normal Lifetimes",
xlab = "Lifetime", ylab = "Density", border = "black")

# Add the theoretical log-normal PDF curve

curve(dlnorm(x, meanlog = mu, sdlog = sigma), col = "red", lwd = 2, add = TRUE)

# Analyze the distribution

summary(random_lifetimes)

# Boxplot for further analysis

boxplot(random_lifetimes, main = "Boxplot of Log-Normal Lifetimes", ylab = "Lifetime",
xlab = "Log-Normal Lifetimes")

5. Exercises:

1. Let us consider a component that follows a Weibull failure rate. Simulate this case and then compare it
to an exponential distribution in R. Analyze the differences in the failure distributions.

Results:

Min. 1st Qu. Median Mean 3rd Qu. Max.

Values 0.122 2.702 4.223 4.451 5.856 13.850

Table 1: Summary Statistics of Weibull Failure Distribution

Min. 1st Qu. Median Mean 3rd Qu. Max.

Values 0.00448 1.47672 3.57407 5.16528 7.08636 32.10942

Table 2: Summary Statistics of Exponential Failure Distribution

Analysis and Interpretation

• The minimal values indicate that the Exponential distribution allows for very early failures compared
to the Weibull distribution.

• The 1st Quartile of the Weibull (2.702) is higher than that of the Exponential (1.47672), indicating
that 25% of Weibull failures occur later than those in the Exponential distribution.

7
 • The median for the Weibull distribution (4.223) is also higher than that for the Exponential distri-
bution (3.57407), suggesting that half of the Weibull failures occur later than half of the Exponential
failures.

• The 3rd Quartile for the Weibull distribution (5.856) is lower than that of the Exponential distri-
bution (7.08636), again reflecting the tendency of Weibull failures to occur earlier than those of the
Exponential distribution.

• The maximum value of the Exponential distribution (32.10942) is significantly higher than that of
the Weibull (13.850), indicating that the Exponential distribution has a wider range of possible
failure times.

• The mean failure time of the Weibull distribution (4.451) is lower than that of the Exponential
distribution (5.16528). This reflects the Weibull’s shape parameter, which leads to fewer high-value
outliers, resulting in earlier failures compared to the Exponential distribution.

Figure 4: Comparison between two failure distributions

# Parameters set up
shape_weibull = 2 # Shape parameter
scale_weibull = 5 # Scale parameter

rate_exponential = 1 / scale_weibull # Rate parameter

# Simulate 1000 random failures from Weibull distribution

[Link](123) # For reproducibility
weibull_failures = rweibull(1000, shape = shape_weibull, scale = scale_weibull)

# Simulate 100 random failures from Exponential distribution

exponential_failures = rexp(1000, rate = rate_exponential)

# Histogram plot of the failure distributions

8
 par(mfrow = c(1, 1))

hist(weibull_failures, breaks = 20, probability = TRUE,

col = "lightblue", main = "Weibull failure distributions",
xlab = "Failure Time", ylab = "Density", border = "black")
curve(dweibull(x, shape = shape_weibull, scale = scale_weibull),
col = "red", lwd = 2, add = TRUE)

# Histogram for Exponential distribution

hist(exponential_failures, breaks = 20, probability = TRUE,
col = "lightgrey", main = "Exponential failure distributions",
xlab = "Failure Time", ylab = "Density", border = "black", ylim = c(0,0.18))
curve(dexp(x, rate = rate_exponential),
col = "red", lwd = 2, add = TRUE)

# Analyze the differences in the distributions

summary_weibull = summary(weibull_failures)
summary_exponential = summary(exponential_failures)

cat("Summary of Weibull Failures:\n")

print(summary_weibull)

cat("\nSummary of Exponential Failures:\n")

print(summary_exponential)

2. Generate random variables from a distribution with a given failure rate and plot the failure times.
Ans: On this note, we have already solved it. For the readers, you may try to generate random variables
from log-normal with different sets of parameter values.

3. Fit a Weibull distribution to a set of failure data and estimate the parameters β and η.
Ans: R code

# Load necessary library

library(MASS)

# Generate a random dataset of length 100 from a Weibull distribution

[Link](123) # For reproducibility
beta = 2 # True shape parameter
eta = 5 # True scale parameter
failure_data = rweibull(100, shape = beta, scale = eta)

# Fit the Weibull distribution to the generated data

fit = fitdistr(failure_data, "weibull", start = list(shape = 1, scale = 1))

# Extract the estimated parameters

beta_est = fit$estimate[’shape’]
eta_est = fit$estimate[’scale’]

# Print the estimated parameters

cat("Estimated Parameters:\n")
cat("Shape (beta):", beta_est, "\n")

9
 cat("Scale (eta):", eta_est, "\n")

# Plot the histogram and the fitted Weibull density function

hist(failure_data, breaks = 10, probability = TRUE,
col = "skyblue", main = "Fitted Weibull Distribution",
xlab = "Failure Time", ylab = "Density", border = "white")

# Add the fitted Weibull density function

curve(dweibull(x, shape = beta_est, scale = eta_est),
col = "red", lwd = 2, add = TRUE)

Figure 5: Fitting of Weibull distribution on failure data

4. Take different types of failure rates for different products and compare their failure rates using exponential,
Weibull, and log-normal distributions.
Ans: Left for the readers.

10