# Routine for Multi-Stress Lognormal Plotting Example
# Coded by Reuel Smith 2015-2017
# ========================================================================
# Example Figure 4.7
# Consider the set of data provided in Table 4.1 that represents the ALT
# data of a component in which temperature is the stress-causing agent of
# failure. Further assume that it has already been shown that the scatter
# in times to failure data listed in Table 4.1 can be best modeled by the
# lognormal distribution, and that the Arrhenius model (exponential) is the
# most appropriate life-stress relationship for this particular ALT.
#
# Temperature # of Tested Units Recorded Failure Times (hrs)
# =========== ================= ======================================
# 150 C 3 2350 2560 2980
# 200 C 9 220 250 330 370 380 460 460 510 610
#
# Solution:
# This routine will generate the Lognormal probability plots for the two
# stress levels 150 C (423 K) and 200 C (472 K)
# ========================================================================
# PART 1: Initialize data and constants
# ========================================================================
# sorting the times to failure in ascending order
# CHANGE these values to process your own data
TTF1 <- sort(c(2350,2560,2980))
TTF2 <- sort(c(220,250,330,370,380,460,460,510,610))
# set the length of the TTF vector or the number of failures
n1 <- length(TTF1)
n2 <- length(TTF2)
# CDF or unreliability for probability plot at Stress Levels 1, 2, and 3
i1 <- c(1:n1) # uncensored data points from 1 to n1
i2 <- c(1:n2) # uncensored data points from 1 to n2
# 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)
F1 <- 100*((i1-0.375)/(n1+0.25))
F2 <- 100*((i2-0.375)/(n2+0.25))
# Lognormal scale conversion of the Failure Percent data
F1B <- qnorm(F1/100,mean=0,sd=1)
F2B <- qnorm(F2/100,mean=0,sd=1)
# reliability for probability plot
R1 <- (100-F1)/100
R2 <- (100-F2)/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.
yfit1 <- F1
yfit2 <- F2
xfit1 <- log(TTF1)
xfit2 <- log(TTF2)
# Calculates the slope and the intercept needed to
# extract the Lognormal parameters
pb2_1 <- lm(yfit1 ~ poly(xfit1, 1, raw=TRUE))
pb2_2 <- lm(yfit2 ~ poly(xfit2, 1, raw=TRUE))
intercept1 <- summary(pb2_1)$coefficients[2,1]
intercept2 <- summary(pb2_2)$coefficients[2,1]
slope1 <- summary(pb2_1)$coefficients[1,1]
slope2 <- summary(pb2_2)$coefficients[1,1]
# the log mean of the data
meant1 <- (50-slope1)/intercept1
meant2 <- (50-slope2)/intercept2
# the 84% value for LN(t)
t84_1 <- (84-slope1)/intercept1
t84_2 <- (84-slope2)/intercept2
# the log-sigma
sigmat1 <- t84_1-meant1
sigmat2 <- t84_2-meant2
# CHANGE NROW if you have a different number of TTF data sets
lognresults <- matrix(c(meant1,sigmat1,meant2,sigmat2), nrow = 2, ncol = 2, byrow = TRUE,dimnames = list(c("Lognormal Parameters 1","Lognormal Parameters 2"),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
yfit2 <- c(qnorm(0.5,mean=0,sd=1),qnorm(0.84,mean=0,sd=1))
xfit2_1 <- c(meant1,meant1+sigmat1)
xfit2_2 <- c(meant2,meant2+sigmat2)
pb2_B1 <- lm(yfit2 ~ poly(xfit2_1, 1, raw=TRUE))
pb2_B2 <- lm(yfit2 ~ poly(xfit2_2, 1, raw=TRUE))
min_mu_sig1 <- summary(pb2_B1)$coefficients[2,1]
min_mu_sig2 <- summary(pb2_B2)$coefficients[2,1]
siginv1 <- summary(pb2_B1)$coefficients[1,1]
siginv2 <- summary(pb2_B2)$coefficients[1,1]
ttfc1 <- c(exp((qnorm(0.001,mean=0,sd=1)-siginv1)/min_mu_sig1),exp((qnorm(0.999,mean=0,sd=1)-siginv1)/min_mu_sig1))
ttfc2 <- c(exp((qnorm(0.001,mean=0,sd=1)-siginv2)/min_mu_sig2),exp((qnorm(0.999,mean=0,sd=1)-siginv2)/min_mu_sig2))
totttfc <- c(ttfc1,ttfc2)
# Computes the upper and lower bound for the TTF axis in terms of log-time
signs1 <- c(floor(log10(min(totttfc))):ceiling(log10(max(totttfc))))
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. Also ADD or SUBTRACT stress
# levels as needed according to the setup
# STRESS LEVEL 1 PLOT
plot(TTF1, F1B, log="x", col="blue",
xlab="Time (hours)", ylab="Percent Failure", pch=16,
xlim=c(min(logtimes1), max(logtimes1)), ylim=c(min(fcB), max(fcB)), axes=FALSE)
lines(ttfc1, fcB, col="blue")
# STRESS LEVEL 2 PLOT
points(TTF2, F2B, col="red", pch=17)
lines(ttfc2, fcB, col="red")
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("Stress Level 1 - 423 deg K", "Stress Level 2 - 473 deg K"),
col=c("blue", "red"), pch=c(16,17), cex=0.8,
text.font=4, bg='white')
lognresults
Top