###############################################
## Author: Joshua M. Tebbs
## Date: 11 December 2011
## Update: 30 December 2017
## STAT 509 course notes: R Code
## Chapter 5
###############################################

# Figure 5.1
# Page 57
# Weibull PDFs
t = seq(0,25,0.01)
# Plot Weibull pdf with beta = 2 and eta = 5
plot(t,dweibull(t,2,5),type="l",lty=1,xlab="t",ylab="PDF")
# Add other pdfs
lines(t,dweibull(y,2,10),lty=4)
lines(t,dweibull(y,3,10),lty=8)
abline(h=0)
# Add legend
legend(15,0.15,lty = c(1,4,8), 
	c(
	expression(paste(beta, "=2, ", eta, "=5")),
	expression(paste(beta, "=2, ", eta, "=10")),
	expression(paste(beta, "=3, ", eta, "=10"))
	))


# Example 5.1
# Page 58-59
# Weibull PDF and CDF
# Plot PDF
t = seq(0,30,0.01)
pdf = dweibull(t,2,10)
plot(t,pdf,type="l",xlab="t",ylab="PDF",cex.lab=1.25)
abline(h=0)
# Plot CDF
cdf = pweibull(t,2,10)
plot(t,cdf,type="l",xlab="t",ylab="CDF",ylim=c(0,1),cex.lab=1.25)
abline(h=0)


# Example 5.2
# Weibull hazard functions
# Page 61-62
t = seq(0,4,0.01)
haz.1 = 3*t^2
haz.2 = 1.5*t^(0.5)
haz.3 = 1+t*0
haz.4 = 0.5*t^(-0.5)
# Plots below were made separately
plot(t,haz.1,type='l',xlab="t",ylab="HAZARD",cex.lab=1.25)
abline(h=0,lty=2)
plot(t,haz.2,type='l',xlab="t",ylab="HAZARD",cex.lab=1.25)
abline(h=0,lty=2)
plot(t,haz.3,type='l',xlab="t",ylab="HAZARD",cex.lab=1.25)
abline(h=0,lty=2)
plot(t,haz.4,type='l',xlab="t",ylab="HAZARD",ylim=c(0,5),cex.lab=1.25)
abline(h=0,lty=2)


# Example 5.3
# Pages 63-67
# Need to install (then load) fitdistrplus package in R
# Enter data
cart.data = c(3.9,4.2,5.4,6.5,7.0,8.8,9.2,11.4,14.3,15.1,15.3,15.5,17.9,18.0,19.0,19.0,23.9,24.8,26.0,34.2)
# Fit Weibull model
fitdist(cart.data,"weibull") 

# Page 65
# Plot PDF, CDF, survivor, and hazard
# I constructed each function separately
# Use maximum likelihood estimates beta.hat, eta.hat (obtained from fitdist output)
beta.hat = 1.99
eta.hat = 16.94
t = seq(0,50,0.01)
# Plot PDF
pdf = dweibull(t,beta.hat,eta.hat)
plot(t,pdf,type="l",xlab="t",ylab="PDF",cex.lab=1.25)
abline(h=0)
# Plot CDF
cdf = pweibull(t,beta.hat,eta.hat)
plot(t,cdf,type="l",xlab="t",ylab="CDF",ylim=c(0,1),cex.lab=1.25)
abline(h=0)
# Plot survivor function
survivor = 1-pweibull(t,beta.hat,eta.hat)
plot(t,survivor,type="l",xlab="t",ylab="SURVIVOR",ylim=c(0,1),cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)
# Plot hazard function
hazard = pdf/survivor
plot(t,hazard,type="l",xlab="t",ylab="HAZARD",cex.lab=1.25)
abline(h=0)

# Page 67
# Weibull qqplot
# Need to install (then load) car package in R 
qqPlot(cart.data,distribution="weibull",shape=beta.hat,scale=eta.hat,xlab="Weibull quantiles",ylab="Cart data",pch=16)