# R & WinBUGS Routine to Perform Bayesian Parameter Estimation
# Coded by Reuel Smith 2015-2017
#
# NOTE:
# Install and/or Load the R2WINBUGS Package prior to using this routine
# (Or R2OPENBUGS depending on what program you use. WinBUGS or OpenBUGS
# must be installed on the computer prior to using this routine.
# ========================================================================
# Example Problem 5.2
# Consider the following example involving a test of a particular metal
# alloy. The sliding test was conducted over a range of different applied
# weights in order to study the effect of weight and to gain a better
# understanding of the wear mechanism. The accelerated degradation test
# data is presented in Table 5.2 and plotted in log-log axes (Figure 5.8).
# Estimate the parameters of power degradation model using the MLE method,
# assuming no measurement errors. Repeat the problem assuming 10#
# multiplicative measurement error.
#
# ========================================================================
# Weight Cycles (hundreds)
# ==================================================================
# Unit 2 5 10 20 50 100 200 500
# ========================================================================
# 10g 1 3.2 4.1 4.5 4.7 5.8 6.8 7.7 9.6
# 2 2.7 3.4 3.8 3.9 5.4 5.7 6.3 8.4
# 3 2.1 2.7 3.1 3.3 4.0 4.6 5.7 6.6
# 4 2.6 3.5 4.0 4.0 5.2 6.1 6.7 8.5
# 50g 5 7.5 7.8 8.2 10.6 12.6 13.3 12.9 14.8
# 6 7.5 8.1 9.8 10.9 14.8 16.1 17.3 20.2
# 7 7.0 8.9 9.4 11.1 12.4 13.5 16.7 17.3
# 8 7.8 8.9 10.0 11.5 13.7 16.2 16.2 21.0
# 100g 9 12.5 15.4 17.2 20.5 24.1 27.0 29.4 37.9
# 10 11.0 13.9 16.1 18.6 22.2 27.8 31.0 36.6
# 11 13.0 15.1 18.6 20.2 23.9 29.7 31.5 39.6
# 12 11.7 13.7 16.7 17.5 22.3 25.3 32.0 38.2
# ========================================================================
# Solution:
# The example problem computes the mean and sigma for all the unit data.
# This example will only compute the parameter distribution of one of the
# units based on a Bayesian parameter analysis of
# log(D(cycles)) = log(sigma1) + sigma2*log(D(cycles))
# This is expected to follow a lognormal distribution.
# ========================================================================
# Initialize data and constants
# ========================================================================
# Definitions of different variables used in this code
# x(1): sigma 1 parameter
# x(2): sigma 2 parameter
# x(3): standard deviation parameter
# Cycle data (in hundreds of cycles)
N <- c(2,5,10,20,50,100,200,500)
Nfull <- c(1:500)
# Degradation measure width of scar (in microns). This degradation measure
# is used as the unit of damage at cycle N in this problem
munit1 <- c(3.2,4.1,4.5,4.7,5.8,6.8,7.7,9.6)
munit2 <- c(2.7,3.4,3.8,3.9,5.4,5.7,6.3,8.4)
munit3 <- c(2.1,2.7,3.1,3.3,4.0,4.6,5.7,6.6)
munit4 <- c(2.6,3.5,4.0,4.0,5.2,6.1,6.7,8.5)
munit5 <- c(7.5,7.8,8.2,10.6,12.6,13.3,12.9,14.8)
munit6 <- c(7.5,8.1,9.8,10.9,14.8,16.1,17.3,20.2)
munit7 <- c(7.0,8.9,9.4,11.1,12.4,13.5,16.7,17.3)
munit8 <- c(7.8,8.9,10.0,11.5,13.7,16.2,16.2,21.0)
munit9 <- c(12.5,15.4,17.2,20.5,24.1,27.0,29.4,37.9)
munit10 <- c(11.0,13.9,16.1,18.6,22.2,27.8,31.0,36.6)
munit11 <- c(13.0,15.1,18.6,20.2,23.9,29.7,31.5,39.6)
munit12 <- c(11.7,13.7,16.7,17.5,22.3,25.3,32.0,38.2)
m <- c(munit1,munit2,munit3,munit4,munit5,munit6,munit7,munit8,munit9,munit10,munit11,munit12)
# number of parameters in likelihood equation: theta1, theta2, sigma (SD)
n <- 3
# ====================================================================
# SETUP FOR SAMPLING
# ====================================================================
nsample <- 100
niter <- 60000
theta1_L <- -100
theta1_H <- 100
theta2_L <- -100
theta2_H <- 100
sigma_L <- 0
sigma_H <- 1000
t <- N
D <- munit1
nt <- length(t)
data <- list("theta1_L","theta1_H","theta2_L","theta2_H","sigma_L","sigma_H","t","D","nt")
parameters <- c("theta1","theta2","sigma")
inits <- list(list(A = 0.00000000001, p = 1, beta = 1))
# Under MODEL.FILE, change the folder name to where your WinBUGS model
# is located. Under BUG.DIRECTORY, change the folder name to where the
# WinBUGS executable is located.
ADT.sim <- bugs(data,inits,parameters,model.file="C:/Users/ReuelS/Documents/ENRE447/Bayesian_Parameter_Model_Example_5_2.odc",n.chains = 1, n.iter = niter, n.burnin=10000, n.thin=1, debug=TRUE, bugs.directory="C:/Users/ReuelS/Desktop/POD Research/WinBUGS/winbugs14/WinBUGS14", DIC=FALSE)
attach.bugs(ADT.sim)
#detach.bugs()
flush.console()
fbest <- exp(mean(theta1) + mean(theta2)*log(Nfull))
# The probability plot
plot(N, munit1, col="blue",
xlab="Fatigue Cycles", ylab="Weight Loss (micrograms)", pch=19,
xlim=c(0, 500), ylim=c(0, 10), axes=FALSE)
lines(Nfull, fbest, lty=4, col="blue")
axis(1, at=NULL, las=2, cex.axis=0.7)
axis(2, at=NULL, las=2, cex.axis=0.7)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = "dotted",lwd = par("lwd"), equilogs = TRUE)
# 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(0, 10, legend=c("Degradation Data","Mean Degradation Curve"),
col=c("blue", "blue"), pch=c(19,-1), lty=c(0,4), cex=0.8,
text.font=4, bg='white')
tableresults <- matrix(c(min(theta1),max(theta1),mean(theta1),sd(theta1),min(theta2),max(theta2),mean(theta2),sd(theta2),min(sigma),max(sigma),mean(sigma),sd(sigma)), nrow = 3, ncol = 4, byrow = TRUE,dimnames = list(c("theta1", "theta2", "sigma"),c("min", "max", "mean", "SD")))
Top