# R-Script to Perform MLE on a Step-Stress Testing Example
# Coded by Taotao Zhou and Reuel Smith 2015-2017
# ========================================================================
# Example Problem 4.17
# For the data given in Table 4.8,
# ========================================================================
# Step Stress Test Time Test Result (time of failure occurrence)
# (volts) (hrs)
# ==== ======= ========= ==========================================
# 1 2 250 0 fails in 11
# 2 3 100 3 fails in 11 (30; 60; 80)
# 3 4 20 3 fails in 8 (2; 10; 16)
# 4 5 10 3 fails in 5 (1; 4; 8)
# 5 6 10 2 fails in 2 (1; 5)
# 6 7 10 0 fails in 0
# ========================================================================
# stresses are
#
# S_1=2, 0<t<=250
# S_2=3, 250<t<=350
# S_3=4, 350<t<=370
# S_4=5, 370<t<=380
# S_5=6, 380<t<=390
# S_6=7, 390<t<=Inf
#
# Using the log-likelihood function in Equation (4.139) estimate the
# parameters beta, K, and n.
# ========================================================================
# Solution:
# The negative log-likelihood loglik is input directly into the R-script
# and minimized (the actual log-likelihood is maximized) to find the MLE
# parameter estimates for beta, K, and n.
# ========================================================================
# Initialize data and constants
# ========================================================================
# Definitions of different variables used in this code
# theta[1]: beta parameter
# theta[2]: K parameter
# theta[3]: n parameter
# Stress Levels 1 through 6 in Volts
V <- c(2,3,4,5,6,7)
# The end time of each of the Stress Levels 1 through 5 (in hours)
t_end <- c(250,350,370,380,390)
# The given ungrouped failure times
ttf <- c(280,310,330,352,360,366,371,374,378,381,385)
# ========================================================================
# The next step is to define the the failure and survival times as well as
# failure and survival stress levels for each of the six step stress
# levels. Note that all instances where the failure time vector (TTF) and
# survival time vector (TTS) is left as 0, are not added to the operation
# because the log-likelihood only counts the times failed and succeded
# through all the steps.
# ========================================================================
# STEP 1: 2V; 250hr test time; 0 fails in 11
# Failure Times in hours (0)
# N/A
# Survival times in hours (11)
TTS1 <- t_end[1]*c(1,1,1,1,1,1,1,1,1,1,1)
# Failure Stress (V)
# N/A
# Survival Stress (V)
VS1 <- V[1]*c(1,1,1,1,1,1,1,1,1,1,1)
# STEP 2: 3V; 100hr test time; 3 fails in 11 (30; 60; 80)
# Failure Times in hours (3)
TTF2 <- c(ttf[1:3])
# Survival times in hours (8)
TTS2 = t_end[2]*c(1,1,1,1,1,1,1,1)
# Failure Stress (V)
VF2 <- V[2]*c(1,1,1)
# Survival Stress (V)
VS2 = V[2]*c(1,1,1,1,1,1,1,1)
# STEP 3: 4V; 20hr test time; 3 fails in 8 (2; 10; 16)
# Failure Times in hours (3)
TTF3 <- c(ttf[4:6])
# Survival times in hours (5)
TTS3 <- t_end[3]*c(1,1,1,1,1)
# Failure Stress (V)
VF3 <- V[3]*c(1,1,1)
# Survival Stress (V)
VS3 <- V[3]*c(1,1,1,1,1)
# STEP 4: 5V; 10hr test time; 3 fails in 5 (1; 4; 8)
# Failure Times in hours (3)
TTF4 <- c(ttf[7:9])
# Survival times in hours (2)
TTS4 <- t_end[4]*c(1,1)
# Failure Stress (V)
VF4 <- V[4]*c(1,1,1)
# Survival Stress (V)
VS4 <- V[4]*c(1,1)
# STEP 5: 6V; 10hr test time; 2 fails in 2 (1; 5)
# Failure Times in hours (2)
TTF5 <- c(ttf[10:11])
# Survival times in hours (8)
# N/A
# Failure Stress (V)
VF5 = V[5]*c(1,1)
# Survival Stress (V)
# N/A
# STEP 6: 7V; 10hr test time; 0 fails in 0
# Failure Times in hours (0)
# N/A
# Survival times in hours (0)
# N/A
# Failure Stress (V)
# N/A
# Survival Stress (V)
# N/A
# ========================================================================
# SETUP FOR TAU (tau1 through tau5 for stress levels 2 through 6
# This establishes the equivalent time (tau_(i-1)) of step i-1 if the item
# is operated at step stress level i. Tau_1 follows Equation 4.135 in the
# textbook while all other tau values follow Equation 4.136.
# ========================================================================
# tau for step 2
tau1 <- function(n) {
t_end[1]*(V[1]/V[2])^n
}
# tau for step 3
tau2 <- function(n) {
(t_end[2]-t_end[1]+tau1(n))*(V[2]/V[3])^n
}
# tau for step 4
tau3 <- function(n) {
(t_end[3]-t_end[2]+tau2(n))*(V[3]/V[4])^n
}
# tau for step 5
tau4 <- function(n) {
(t_end[4]-t_end[3]+tau3(n))*(V[4]/V[5])^n
}
# tau for step 6
tau5 <- function(n) {
(t_end[5]-t_end[4]+tau4(n))*(V[5]/V[6])^n
}
# ========================================================================
# SETUP FOR ADJUSTED TIME VECTORS
# The procedure for setting up the adjusted failure time vector ti, and the
# adjusted survival time vector tj follows the expression given in Equation
# 4.134. The vector tvect is the t_(i-1) value, the end time of step i-1,
# and tau is the tau_(i-1) value. Each failure time and survival time is
# assigned a t_(i-1) and tau_(i-1) according to the step.
tvecti <- c(t_end[1]*c(1,1,1),t_end[2]*c(1,1,1),t_end[3]*c(1,1,1),t_end[4]*c(1,1))
tvectj <- c(c(0,0,0,0,0,0,0,0,0,0,0),t_end[1]*c(1,1,1,1,1,1,1,1),t_end[2]*c(1,1,1,1,1),t_end[3]*c(1,1))
taui <- function(n) {
c(tau1(n),tau1(n),tau1(n),tau2(n),tau2(n),tau2(n),tau3(n),tau3(n),tau3(n),tau4(n),tau4(n))
}
tauj <- function(n) {
c(c(0,0,0,0,0,0,0,0,0,0,0),tau1(n),tau1(n),tau1(n),tau1(n),tau1(n),tau1(n),tau1(n),tau1(n),tau2(n),tau2(n),tau2(n),tau2(n),tau2(n),tau3(n),tau3(n))
}
ti <- function(n) {
c(TTF2,TTF3,TTF4,TTF5) - tvecti + taui(n)
}
tj <- function(n) {
c(TTS1,TTS2,TTS3,TTS4) - tvectj + tauj(n)
}
# The failure and survival stress levels are assigned a vector of the same
# sizes as the adjusted time vectors for the MLE operation.
Vi <- c(VF2,VF3,VF4,VF5)
Vj <- c(VS1,VS2,VS3,VS4)
# SETUP EQUATION AND SOLVE MLE
# The equation for MLE analysis, the log-likelihood needs to be flipped to
# the negative since the NLM function is a minimization function. The
# initial estimate c(x1,x2,x3) should be checked in advance to see if it is
# a viable estimate or as close zero as possible. A reasonable first guess
# may be acquired by using the Plotting Step-Stress method defined in
# section 4.16.1. Randomly estimating c(x1,x2,x3) will almost always
# result in a poor estimate.
loglik <- function(theta) {
-(sum(log(theta[1]*theta[2]*(Vi^theta[3])*(theta[2]*(Vi^theta[3])*ti(theta[3]))^(theta[1]-1))) - sum((theta[2]*(Vi^theta[3])*ti(theta[3]))^theta[1]) - sum((theta[2]*(Vj^theta[3])*tj(theta[3]))^theta[1]))
}
nlm(loglik, theta <- c(2.5996,0.00015462,1.8308), hessian=TRUE)
Top