Weibull Analysis with WeibullR

Paul Govan

The Weibull Distribution

The Weibull distribution is the most common distribution in the field of life data analysis

Why? Because the Weibull is a flexible distribution that can fit many different types of data

The Normal Distribution

Take the following normal distribution

x <- seq(0, 10, by = .1)
y <- dnorm(x, mean = 5, sd = 1.5)
plot(x,y, type = "l")

The Weibull vs Normal

The Weibull can mimic the normal and many others

y <- dweibull(x, shape = 3.678, scale = 5.54)
plot(x,y, type = "l")

The Cumulative Distribution Function

R(t) gives the cumulative probability of survival prior to time t while F(t) gives the cumulative probability of failure prior to time t

\[ R(t) = 1 - F(t) = e^-(t/\eta)^\beta \]

Where \(\beta\) is the shape parameter. In life analysis \(\beta\) has a physical interpretation:

\(\beta\) < 1 - the failure rate is decreasing (may indicate infant mortality)
\(\beta\) = 1 - the failure rate is constant (random failures or possible mix of failure modes)
\(\beta\) > 1 - the failure rate is increasing (may indicate wear out)

And \(\eta\) is the scale parameter or the quantile at 63.2% probability

The Weibull Parameters

Let’s vary the parameters of the Weibull distribution to see these relationships graphically

Varying Beta

Varying Beta changes the shape of the cumulative distribution function

x <- seq(0, 50, by = .1)
plot(pweibull(x, shape = 0.5, scale = 10), col = "blue", type = "l", xlab = "x", ylab = "y") 
lines(pweibull(x, shape = 1, scale = 10), col = "red") 
lines(pweibull(x, shape = 2, scale = 10), col = "green") 
legend(300, 0.6, legend=c("shape = 0.5, scale = 10", "shape = 1, scale = 10", "shape = 2, scale = 10"),col=c("blue", "red", "green"), lty=1, cex=1.2)

Varying Eta

Varying Eta moves the weibull distribution along the x-axis

x <- seq(0, 50, by = .1)
plot(pweibull(x, shape = 1, scale = 5), col = "blue", type = "l", xlab = "x", ylab = "y") 
lines(pweibull(x, shape = 1, scale = 10), col = "red") 
lines(pweibull(x, shape = 1, scale = 20), col = "green") 
legend(300, 0.6, legend=c("shape = 1, scale = 5", "shape = 1, scale = 10", "shape = 1, scale = 20"), col=c("blue", "red", "green"), lty=1, cex=1.2)

WeibullR

Getting Started with WeibullR

For this tutorial, we will use WeibullR: an R package for life data analysis in the tradition of Walodi Weibull

First, check if WeibullR is installed in R and install if not

if(!(require(WeibullR))){install.packages('WeibullR')}

Let’s run a basic example using the quick fit method

failures<-c(30, 49, 82, 90, 96)
fit<-MLEw2p(failures)
knitr::kable(fit, "html", col.names = "param")

	param
Eta	77.811667
Beta	3.201164
LL	-23.181594

Data Censoring

Right Censored Data

Right censored data includes suspensions or units that have operated for a period of time without failure

Let’s add right censored data to the previous example and build a basic probability plot

suspensions<-c(100, 45, 10)
fit<-MLEw2p(failures, suspensions, bounds=TRUE, show=TRUE)

Interval Censored Data

Interval censored data is most commonly in the form of inspection data, where the failure occurred between two inspection times, the latest known working time (left entry), and the inspection time when the failure was detected (right entry)

To create an interval censored model, let’s use the inspection data from Silkworth, 2020

inspection_data<-data.frame(left=c(0, 6.12, 19.92, 29.64, 35.4, 39.72, 45.32, 52.32),
                            right=c(6.12, 19.92, 29.64, 35.4, 39.72, 45.32, 52.32, 63.48),
                            qty=c(5, 16, 12, 18, 18, 2, 6, 17))

Add suspension data for units surviving until the last inspection date

suspensions<-data.frame(time=63.48, event=0, qty=73)

Interval Censored Model

Add a fit and plot the results

obj1<-wblr(suspensions, interval = inspection_data)
obj1<-wblr.fit(obj1, method.fit="mle", col="red") 
obj1<-wblr.conf(obj1, method.conf="fm", lty=2)
plot(obj1)

Grouped Data

Life data is often compiled into groups (i.e. a group of units failed or survived at a certain time) Notice how the previous example contained groups in the qty column

suspensions<-data.frame(time=63.48, event=0, qty=73)
knitr::kable(suspensions, "html")

time	event	qty
63.48	0	73

In this case, there are 73 suspensions at time 63.48

Parameter Estimation Methods

Maximum Likelihood vs Rank Regression

Median Rank Regression (MRR) is a least squares analysis based on minimizing the squared distance from the cumulative probability to the fitted model
Maximum Likelihood Estimation (MLE) is based on estimating the parameters of the distribution by maximizing the likelihood function of the observed data

WeibullR includes functions for both methods, such as the MRRw2p, MRRw3p functions for MRR and the MLEw2p, MLEw3p functions for MLE

Other Weibull Models

The WeiBayes Model

A WeiBayes or 1P Weibull has a fixed \(\beta\) or shape parameter based on prior knowledge or experience

Weibayes models are typically used when failure data is limited

A WeiBayes Example

Let’s rerun the previous example as a WeiBayes model with a Beta of 2

y <- dnorm(x, mean = 5, sd = 1)
suspensions<-c(100, 45, 10)
eta <- weibayes(failures, suspensions, beta=2)
x <- seq(0, 200, by = .1)
plot(pweibull(x, shape = 2, scale = eta), col = "blue", type = "l", xlab = "Time to Failure", ylab = "Unreliability (%)")

The 3-Parameter Weibull Model

The 3P Weibull has an additional parameter t₀ to represent a failure free period (i.e. a component cannot fail until after time t₀)

\[ R(t) = 1 - F(t) = e^-((t-t~0~)/\eta)^\beta \]

A 3P Weibull Example

Let’s rerun the previous example as a 3P Weibull

failures<-c(30, 49, 82, 90, 96)
suspensions<-c(100, 45, 10)
fit<-MLEw3p(failures, suspensions, bounds=TRUE, show=TRUE)

Competing Failure Modes

A system or component may have multiple different failure modes (e.g. corrosion, fatigue, overload)

Each failure mode has a different distribution of time to failure

A basic Weibull with competing failure modes tends to -

Not fit the data well
Look random (Beta of 1)

Competing FM Example

First, create some Weibull data representing 3 different failure modes

data <- data.frame(
  time = c(rweibull(5, 0.5, 10), rweibull(10, 1, 10), rweibull(5, 2, 10), rweibull(100, 2, 10)),
  event = c(rep(1, 20), rep(0, 100)),
  failure_mode = c(rep("A", 5), rep("B", 10), rep("C", 5), rep("", 100))
  )

Competing FM Example (Cont’d)

Fit an overall Weibull

obj <- wblr.conf(wblr.fit(wblr(data)))
plot(obj, col="darkgreen")

Competing FM Example (Cont’d)

Now create a separate data set for each failure mode

dat1 <- data
dat2 <- data
dat3 <- data
dat1$event[dat1$failure_mode != "A"] <- 0
dat2$event[dat2$failure_mode != "B"] <- 0
dat3$event[dat3$failure_mode != "C"] <- 0

Competing FM Example (Cont’d)

Finally, fit a wblr object for each failure mode and plot the results

obj1 <- wblr.fit(wblr(dat1, col="blue"))
obj2 <- wblr.fit(wblr(dat2, col="red"))
obj3 <- wblr.fit(wblr(dat3, col="orange"))
plot.wblr(list(obj1, obj2, obj3), is.plot.legend = FALSE)

Multi-Plots

To build a plot with multiple Weibull models, start by creating a wblr object

failures<-c(30, 49, 82, 90, 96)
obj1<-wblr(failures)

Add a default fit that appears as a red line on the plot

obj1<-wblr.fit(obj1, col="red")

Add default confidence bounds using a single width line

obj1<-wblr.conf(obj1, lwd=1)

Now create a second wblr object

failures<-c(30, 49, 82, 90, 96)
suspensions<-c(100, 45, 10)
obj2<-wblr.conf(wblr.fit(wblr(failures, suspensions, col="purple"),),lwd=1)

Multi-Plots (Cont’d)

Add the 2 wblr objects to a list and plot both objects in a single chart

plot.wblr(list(obj1, obj2))

Wrapping Up

To Get Help

To get help with a specific function, type a question mark before the function name

?MRRw2p

For more help with the WeibullR package

help(package="WeibullR")

References

Silkworth D, Symynck J (2022). WeibullR: Weibull Analysis for Reliability Engineering. R package version 1.2.1, https://CRAN.R-project.org/package=WeibullR.
Silkworth, David. (2020). WeibullR: An R Package for Weibull Analysis for Reliability Engineers. 43-53. https://doi.org/10.35566/isdsa2019c3.