# Lognormal 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 Lognormal solution for the
# parameters log mean and log sigma, and the Lognormal 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)
# CDF or unreliability for probability plot at Stress Levels 1, 2, and 3
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:
# MIDPOINT - (i - 0.5)/n
# MEAN - i/(n + 1)
# MEDIAN - (i - 0.3)/(n + 0.4)
F <- 100*((i-0.375)/(n+0.25))
# Lognormal scale conversion of the Failure Percent data
FB <- qnorm(F/100,mean=0,sd=1)
# 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
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 Lognormal plot scale
# Upper and lower bounds of the Percent Failure axis in percent
fc <- c(.1,99.9)
# Lognormal scale conversion of the upper and lower bounds of the Percent
# Failure axis
fcB <- qnorm(fc/100,mean=0,sd=1)
# 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)
# Lognormal scale conversion of Percent Failure ticks cited between the
# lower and upper bounds
Pticks <- qnorm(Pticks1/100,mean=0,sd=1)
Pticks1label <- c(0.1,0.2,"","",0.5,"","","","",1,2,"","",5,"","","","",10*c(1:9),95,99,99.9)
# ========================================================================
# PART 2: Solve for Lognormal Parameters log mean and log sigma
# ========================================================================
# Use a polynomial fit (polyfit) to draw the best fit line (fc) through the
# known times to failure TTF from F(t) = 1% to F(t) = 99%
# polynomial fit for LN(TTF) and F
#
# 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 <- F
xfit <- log(TTF)
# Calculates the slope and the intercept needed to
# extract the Lognormal parameters
pb2 <- lm(yfit ~ poly(xfit, 1, raw=TRUE))
intercept <- summary(pb2)$coefficients[2,1]
slope <- summary(pb2)$coefficients[1,1]
# the log mean of the data
meant <- (50-slope)/intercept
# the 84% value for LN(t)
t84 <- (84-slope)/intercept
# the log-sigma
sigmat <- t84-meant
lognresults <- matrix(c(meant,sigmat), nrow = 1, ncol = 2, byrow = TRUE,dimnames = list(c("Lognormal Parameters"),c("mean_t", "sigma_t")))
# ========================================================================
# PART 3: Lognormal probability plot
# ========================================================================
# Plot points at 0.1% and 99.9%
# These have a different polynomial pair than that which was used to find
# the meant and sigmat
xfit2 <- c(meant,meant+sigmat)
yfit2 <- c(qnorm(0.5,mean=0,sd=1),qnorm(0.84,mean=0,sd=1))
pb2_B <- lm(yfit2 ~ poly(xfit2, 1, raw=TRUE))
min_mu_sig <- summary(pb2_B)$coefficients[2,1]
siginv <- summary(pb2_B)$coefficients[1,1]
ttfc <- c(exp((qnorm(0.001,mean=0,sd=1)-siginv)/min_mu_sig),exp((qnorm(0.999,mean=0,sd=1)-siginv)/min_mu_sig))
# 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
# 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[2]), 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 the (X,Y) 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], qnorm(0.99,mean=0,sd=1), legend=c("Data", "Lognormal best fit line"),
col=c("blue", "blue"), pch=c(16,-1), lty=c(0,1), cex=0.8,
text.font=4, bg='white')
lognresults