Lognormal Plotting of Single Data (Example Problem 5.1), R

# 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

# 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 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

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 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

Top