# Weibull Plotting (Single) Example
# Coded by Reuel Smith 2015-2017
# ========================================================================
# Example Problem 5.1
# Five specimens of a new corrosion-resistant material are tested for 240
# hours in a highly corrosive environment. The density of the material is
# 7.6 g/cm3, and the exposed surface area of each specimen is 4.3 cm2. At
# the end of the test period, the measured weight losses in mg were 11.1,
# 10.4, 12.1, 11.4, and 9.8 (assume no measurement error). If a
# degradation of 1 mm or more results in a structural failure, predict the
# failure times for the five specimens. Based on these failure times,
# determine the probability distribution that best represents the life of
# the material.
#
# Solution:
# TTF = 70659, 75415, 64820, 68800, and 80033
# The proposed probability solutions for this problem are Weibull and
# Lognormal. This routine will process the Weibull solution for the
# parameters alpha and beta, and the Weibull fit probability plot
# ========================================================================
# PART 1: Initialize data and constants
# ========================================================================
# sorting the times to failure in ascending order
# CHANGE THESE VALUES to process your own data
TTF <- sort(c(70659,75415,64820,68800,80033))
# set the length of the TTF vector or the number of failures
n <- length(TTF)
# uncensored data points from 1 to n
i <- c(1:n) # uncensored data points from 1 to n
# CDF or Failure Percent for probability plot under the KIMBALL probability
# plotting position: (i - 0.375)/(n + 0.25). Other PLOTTING POSITIONS may
# be assumed as given:
# MIDPOINT - (i - 0.5)/n
# MEAN - i/(n + 1)
# MEDIAN - (i - 0.3)/(n + 0.4)
F <- 100*((i-0.375)/(n+0.25))
# Weibull scale conversion of the Failure Percent data
FB <- log(log(1/(1-F/100)))
# reliability for probability plot
R <- (100-F)/100
# Form a matrix to categorize the results of the above analysis.
# Column 1 displays the TTFs, Column 2 displays the calculated Failure
# Probability, and Column 3 displays the Reliability
# NOTE: Remember that when your data count changes you will have to make additional columns manually.
# Add to the "dimnames" list.
TTFtable <- matrix(c(1:(3*n)), nrow = n, ncol = 3, byrow = TRUE,dimnames = list(c("1", "2", "3", "4", "5"),c("TTF", "Failure Probaility (%)", "Reliability (%)")))
TTFtable[,1] <- TTF
TTFtable[,2] <- F
TTFtable[,3] <- R*100
# Scale conversion to Weibull plot scale
# Upper and lower bounds of the Percent Failure axis in percent
fc <- c(.1,99.9)
# Weibull scale conversion of the upper and lower bounds of the Percent
# Failure axis
fcB <- log(log(1/(1-fc/100)))
# Percent Failure ticks cited between the lower and upper bounds
Pticks1 <- c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,c(1:10),10*c(2:9),95,99,99.9)
# Weibull scale conversion of Percent Failure ticks cited between the lower
# and upper bounds
Pticks <- log(log(1/(1-Pticks1/100)))
Pticks1label <- c(0.1,0.2,0.3,"",0.5,"","","","",1,2,3,"",5,"","","","",10*c(1:9),95,99,99.9)
# ========================================================================
# PART 2: Solve for Weibull Parameters alpha and beta
# ========================================================================
# The Weibull parameters alpha and beta may be computed by assuming a
# relationship for the reliability function R(t) = exp[-(t/alpha)^beta] or
# 1/R(t) = exp[(t/alpha)^beta]. The relation translates to:
#
# ln[ln(1/R(t))] = beta*ln(t) - beta*ln(alpha)
#
# The LM function is used to fit this linear relationship and extract
# the Weibull parameters alpha and beta. See Reliability Engineering and
# Risk Analysis: A Practical Guide Section 3.3.2.2. for method.
yfit <- log(log(1/R))
xfit <- log(TTF)
# Calculates the slope beta and the intercept beta*ln(alpha) needed to
# extract the Weibull parameters
wblparm <- lm(yfit ~ poly(xfit, 1, raw=TRUE))
beta <- summary(wblparm)$coefficients[2,1]
alpha = exp(-summary(wblparm)$coefficients[1,1]/beta)
wblresults <- matrix(c(alpha,beta), nrow = 1, ncol = 2, byrow = TRUE,dimnames = list(c("Wbl Parameters"),c("alpha", "beta")))
# ========================================================================
# PART 3: Weibull probability plot
# ========================================================================
# Plot points at 0.1% and 99.9%
# These have a different polynomial pair than that which was used to find
# the alpha and beta
intercept <- summary(wblparm)$coefficients[1,1]
ttfc <- c(exp((log(log(1./(1-0.001)))-intercept)/beta),exp((log(log(1./(1-0.999)))-intercept)/beta))
# Computes the upper and lower bound for the TTF axis in terms of log-time
signs1 <- c(floor(log10(ttfc[1])):ceiling(log10(ttfc[2])))
logtimes1 <- 10^signs1
Pticks1X <- c(1:(9*length(logtimes1)-8))
Pticks1X[1] <- logtimes1[1]
Pticks1Xlabel <- Pticks1X
# Calculates the tick values used for the final plot
for(i2 in 1:(length(signs1)-1)){
Pticks1X[(9*i2-7):(9*(i2+1)-8)] <- logtimes1[i2]*c(2:10)
Pticks1Xlabel[(9*i2-7):(9*(i2+1)-8)] <- c("","","","","","","","",logtimes1[i2+1])
}
# The probability plot
# DEFAULT XLAB HEADING is in "Time (hours)". When using a different unit of
# time or measurement, CHANGE XLAB HEADING BELOW
plot(TTF, FB, log="x", col="blue",
xlab="Time (hours)", ylab="Percent Failure", pch=16,
xlim=c(10^signs1[1], 10^signs1[length(signs1)]), ylim=c(min(fcB), max(fcB)), axes=FALSE)
lines(ttfc, fcB, col="blue")
axis(1, at=Pticks1X,labels=Pticks1Xlabel, las=2, cex.axis=0.7)
axis(2, at=Pticks,labels=Pticks1label, las=2, cex.axis=0.7)
#Add horizontal grid
abline(h = Pticks, lty = 2, col = "grey")
#Add vertical grid
abline(v = Pticks1X, lty = 2, col = "grey")
# The legend heading
# The first two places signify where the legend will be located. CHANGE the
# legend labels in the LEGEND="" portion. PCH changes the label
# types and LTY changes the line types.
legend(10^signs1[1], log(log(1/(1-.99))), legend=c("Data", "Weibull best fit line"),
col=c("blue", "blue"), pch=c(16,-1), lty=c(0,1), cex=0.8,
text.font=4, bg='white')
wblresults
Top