###############################################
## Author: Joshua M. Tebbs
## Date: 4 August 2025
## Update: 9 August 2025
## STAT 509 course notes: R Code Chapter 5
###############################################

# Use the following command as necessary
par(mar=c(4,4.5,4,3.5)) # adjust plotting region to avoid cutting off vertical axis label 

# Install and load fitdistrplus R package
# See https://cran.r-project.org/web/packages/fitdistrplus/fitdistrplus.pdf for documentation

# Install and load car R package
# See https://cran.r-project.org/web/packages/car/car.pdf for documentation

# Figure 5.1
# Page 74
t = seq(0,25,0.01)
plot(t,dweibull(t,2,5),type="l",lty=1,xlab="t",ylab=expression(f[T](t)),cex.lab=1.25)
lines(t,dweibull(t,2,10),lty=4)
lines(t,dweibull(t,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, "=5")),
	expression(paste(beta, "=3, ", eta, "=10"))
	))


# Figure 5.2
# Page 76
t = seq(0,12000,1)
# Left (PDF)
pdf = dweibull(t,1.5,3000)
plot(t,pdf,type="l",lty=1,xlab="t",ylab=expression(f[T](t)),cex.lab=1.25)
abline(h=0)
# Right (CDF)
cdf = pweibull(t,1.5,3000)
plot(t,cdf,type="l",lty=1,xlab="t",ylab=expression(F[T](t)),cex.lab=1.25)
abline(h=0)


# Figure 5.3
# Page 78
t = seq(0,4,0.01)
plot(t,2.2*t^(1.2),type="l",lty=1,xlab="t",ylab=expression(h[T](t)),ylim=c(0,7.5),cex.lab=1.25)
lines(t,1.5*t^(0.5),lty=4)
lines(t,0.5*t^(-0.5),lty=8)
abline(h=0)
abline(v=0)


# Figure 5.4
# Page 80
t = seq(0,12000,1)
pdf = dweibull(t,1.5,3000)
cdf = pweibull(t,1.5,3000)
plot(t,pdf/(1-cdf),type="l",lty=1,xlab="t",ylab=expression(h[T](t)),cex.lab=1.25)
abline(h=0)


# Example 5.2 (data analysis)
# Page 82
library(fitdistrplus) 
distance.to.failure = c(6700,6950,7820,9120,9660,9820,11310,11690,11850,11880,12140,12200,12870,
		13150,13330,13470,14040,14300,17520,17540,17890,18450,18960,18980,19410,20100,
		20100,20150,20320,20900,22700,23490,26510,27410,27490,27890,28100,30050)
options(digits=3)
fitdist(distance.to.failure,distr="weibull",method="mle")

#################################################################################################
# The code is used to construct Figures 5.5-5.7 and 5.9 below
library(fitdistrplus)
distance.to.failure = c(6700,6950,7820,9120,9660,9820,11310,11690,11850,11880,12140,12200,12870,
		13150,13330,13470,14040,14300,17520,17540,17890,18450,18960,18980,19410,20100,
		20100,20150,20320,20900,22700,23490,26510,27410,27490,27890,28100,30050)
mles = fitdist(distance.to.failure,distr="weibull",method="mle") 
# Extract beta.hat and eta.hat from the fitdist output to construct figures below
beta.hat = mles$estimate[1] 
eta.hat = mles$estimate[2] 
#################################################################################################

# Figure 5.5
# Page 82
t = seq(0,50000,1)
pdf = dweibull(t,beta.hat,eta.hat) 
plot(t,pdf,type="l",lty=1,xlab="t",ylab=expression(f[T](t)),cex.lab=1.25)
abline(h=0)


# Figure 5.6
# Page 83
t = seq(0,50000,1)
# Left (CDF)
cdf = pweibull(t,beta.hat,eta.hat) 
plot(t,cdf,type="l",lty=1,xlab="t",ylab=expression(F[T](t)),cex.lab=1.25)
abline(h=0)
# Right (Survivor)
plot(t,1-cdf,type="l",lty=1,xlab="t",ylab=expression(S[T](t)),cex.lab=1.25)
abline(h=0)


# Figure 5.7
# Page 84
t = seq(0,50000,1)
plot(t,pdf/(1-cdf),type="l",lty=1,xlab="t",ylab=expression(h[T](t)),cex.lab=1.25)
abline(h=0)


# Figure 5.8
# Page 85
distance.to.failure = c(6700,6950,7820,9120,9660,9820,11310,11690,11850,11880,12140,12200,12870,
		13150,13330,13470,14040,14300,17520,17540,17890,18450,18960,18980,19410,20100,
		20100,20150,20320,20900,22700,23490,26510,27410,27490,27890,28100,30050)
bins = seq(0,35000,5000)
hist(distance.to.failure,breaks=bins,xlab="Distance to failure (in km)",ylab="Count",main="",col="lightblue")


# Figure 5.9
# Page 86
library(car)
distance.to.failure = c(6700,6950,7820,9120,9660,9820,11310,11690,11850,11880,12140,12200,12870,
		13150,13330,13470,14040,14300,17520,17540,17890,18450,18960,18980,19410,20100,
		20100,20150,20320,20900,22700,23490,26510,27410,27490,27890,28100,30050)
qqPlot(distance.to.failure,distribution="weibull",shape=beta.hat,scale=eta.hat,
	xlab="Weibull quantiles",ylab="Distance to failure (in km)",pch=16,
	envelope=list(border=TRUE,style="lines"),id=FALSE)


# Figure 5.10
# Page 87
data = rweibull(50,shape=2,scale=200)
qqPlot(data,distribution="weibull",shape=2,scale=200,
	xlab="Weibull quantiles",ylab="Observed data",pch=16,
	envelope=list(border=TRUE,style="lines"),id=FALSE)