###############################################
## Author: Joshua M. Tebbs
## Date: 11 November 2018
## Update: 30 November 2018
## STAT 512 course notes: R Code
## Chapter 6
###############################################

# Figure 6.1
# Page 3
# Left
y = seq(0,1,0.01)
pdf = rep(1,length(y))
plot(y,pdf,type="l",xlab="y",ylab="PDF",xaxt='n',yaxt='n',xlim=c(0,1.5),ylim=c(0,1.5),cex.lab=1.25)
axis(side=1,at=c(0,1),labels=c("0","1"))
axis(side=2,at=c(0,1),labels=c("0","1"))
lines(c(1,1.5),c(0,0),lty=1)
# Right
y = seq(0.001,1,0.001)
h = -log(y)
plot(y,h,type="l",xlab="y",ylab="h(y)",xlim=c(0,1),ylim=c(0,6),xaxt='n',cex.lab=1.25)
axis(side=1,at=c(0,0.5,1),labels=c("0","0.5","1"))
abline(h=0)


# Figure 6.2
# Page 5
# Left
y = seq(0,6,0.01)
pdf = dexp(y,1)
plot(y,pdf,type="l",xlab="y",ylab="PDF",xaxt='n',yaxt='n',xlim=c(0,6),ylim=c(0,0.9),cex.lab=1.25)
axis(side=1,at=c(0),labels=c("0"))
axis(side=2,at=c(0),labels=c("0"))
abline(h=0)
abline(v=0,lty=2)
# Right
y = seq(0,4,0.001)
m = 3
h = y^(1/m)
plot(y,h,type="l",xlab="y",ylab="h(y)",xlim=c(0,4),ylim=c(0,1.8),xaxt='n',yaxt='n',cex.lab=1.25)
axis(side=1,at=c(0),labels=c("0"))
axis(side=2,at=c(0),labels=c("0"))
abline(h=0)


# Figure 6.3
# Page 6
y = seq(0,20,0.01)
plot(y,dweibull(y,shape=1.5,scale=5^(1/1.5)),type="l",lty=1,xlab="u",ylab="PDF",cex.lab=1.25)
lines(y,dweibull(y,shape=2,scale=50^(1/2)),lty=4)
lines(y,dweibull(y,shape=2.5,scale=100^(1/2.5)),lty=8)
abline(h=0)
legend(10,0.20,lty=c(1,4,8),c(expression(paste(m,"=1.5, ",alpha,"=5")),expression(paste(m,"=2.0, ",alpha,"=50")),expression(paste(m,"=2.5, ",alpha,"=100"))))


# Figure 6.4
# Page 7
# Left
y = seq(-3.5,3.5,0.01)
pdf = dnorm(y,0,1)
plot(y,pdf,type="l",xlab="y",ylab="PDF",xaxp=c(-3,3,6),cex.lab=1.25)
abline(h=0)
# Right
u = seq(0,9,0.001)
pdf = dchisq(u,1)
plot(u,pdf,type="l",xlab="u",ylab="PDF",xlim=c(0,9),ylim=c(0,2),cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)


# Figure 6.5
# Page 8
x <- seq(0,2,0.001)
y1 <- rep(2,length(x))
y2 <- rep(0,length(x))
plot(x,y1,type="l",bty="n",xlim=c(0,2),ylim=c(0,2),xaxt='n',yaxt='n',xlab=expression(y[1]),ylab=expression(y[2]),cex.lab=1.25)
polygon(c(x,rev(x)),c(y2,rev(y1)),col="grey")
axis(side=1,at=c(0,2),labels=c("0",""))
axis(side=2,at=c(0,2),labels=c("0",""))


# Figure 6.6
# Page 9
x <- seq(0,3,0.01)
y1 <- 3-x
y2 <- rep(0,length(x))
plot(x,y1,type="l",bty="n",xlim=c(0,4),ylim=c(0,4),xaxt='n',yaxt='n',xlab=expression(y[1]),ylab=expression(y[2]),cex.lab=1.25)
polygon(c(x,rev(x)),c(y2,rev(y1)),col="grey")
axis(side=1,at=c(0,3,5),labels=c("0","u",""))
axis(side=2,at=c(0,3,5),labels=c("0","u",""))


# Figure 6.7
# Page 12
# Left
y = seq(0,1,0.001)
pdf = dbeta(y,2,6)
plot(y,pdf,type="l",xlab="y",ylab="PDF",cex.lab=1.25)
abline(h=0)
# Right
u = seq(0,1,0.001)
pdf = dbeta(u,6,2)
plot(u,pdf,type="l",xlab="u",ylab="PDF",cex.lab=1.25)
abline(h=0)


# Figure 6.8
# Page 13
# Left
y = seq(-pi/2,pi/2,0.01)
pdf = (1+y-y)/pi
matplot(y,pdf,xlab="y",ylab="PDF",type="l",xaxt="n",yaxt="n",xlim=c(-2.5,2.5),ylim=c(0,0.75),cex.lab=1.25)
axis(1, at = c(-pi/2, -pi/4, 0, pi/4, pi/2),labels = expression(-pi/2, -pi/4, 0, pi/4, pi/2))
axis(2, at = c(0, 1/pi),labels = expression(0, 1/pi))
lines(c(-2.5,-pi/2),c(0,0),lty=1)
lines(c(pi/2,2.5),c(0,0),lty=1)
# Right
y = seq(-pi/2,pi/2,0.01)
u = tan(y)
matplot(y,u,xlab="y",ylab="h(y)",type="l",xaxt="n",yaxt="n",xlim=c(-2.5,2.5),ylim=c(-10,10),cex.lab=1.25)
abline(v=0,lty=2)
abline(v=-pi/2,lty=2)
abline(v=pi/2,lty=2)
axis(1, at = c(-pi/2, -pi/4, 0, pi/4, pi/2),labels = expression(-pi/2, -pi/4, 0, pi/4, pi/2))
axis(2, at = c(-10,-5,0,5,10),labels = c("-10","-5","0","5","10"))


# Figure 6.9
# Page 14
u = seq(-8,8,0.001)
pdf = (1/pi)*(1/(1+u^2))
plot(u,pdf,type="l",xlab="u",ylab="PDF",xaxp=c(-8,8,8),cex.lab=1.25)
abline(h=0)


# Figure 6.10
# Page 16
u = seq(0,30,0.01)
plot(u,dlnorm(u,1,1),type="l",lty=1,xlab="u",ylab="PDF",cex.lab=1.25)
lines(u,dlnorm(u,1.5,0.5),lty=4)
lines(u,dlnorm(u,2.5,0.25),lty=8)
abline(h=0)
legend(15,0.20,lty=c(1,4,8),c(expression(paste(mu,"=1, ",sigma,"=1")),
expression(paste(mu,"=1.5, ",sigma,"=0.5")),expression(paste(mu,"=2.5, ",sigma,"=0.25 "))))


# Figure 6.11
# Page 26
# Left
x <- seq(0,2,0.001)
y1 <- rep(2,length(x))
y2 <- rep(0,length(x))
plot(x,y1,type="l",bty="n",xlim=c(0,2),ylim=c(0,2),xaxt='n',yaxt='n',xlab=expression(y[1]),ylab=expression(y[2]),cex.lab=1.25)
polygon(c(x,rev(x)),c(y2,rev(y1)),col="grey")
axis(side=1,at=c(0,2),labels=c("0",""))
axis(side=2,at=c(0,2),labels=c("0",""))
# Right
x <- seq(0,2,0.001)
y1 <- rep(1/4,length(x))
y2 <- rep(0,length(x))
plot(x,y1,type="l",bty="n",xlim=c(0,2),ylim=c(0,2),xaxt='n',yaxt='n',xlab=expression(u[1]),ylab=expression(u[2]),cex.lab=1.25)
polygon(c(x,rev(x)),c(y2,rev(y1)),col="grey")
axis(side=1,at=c(0,2),labels=c("0",""))
axis(side=2,at=c(0,1/4,2),labels=c("0","1",""))


# Figure 6.12
# Page 29
u = seq(-6,6,0.001)
pdf = 0.5*exp(-abs(u))
plot(u,pdf,type="l",xlab=expression(u[2]),ylab="PDF",ylim=c(0,0.5),cex.lab=1.25)
abline(h=0)


# Figure 6.13
# Page 30
# Left
x <- seq(0,2,0.001)
y1 <- rep(2,length(x))
y2 <- rep(0,length(x))
plot(x,y1,type="l",bty="n",xlim=c(0,2),ylim=c(0,2),xaxt='n',yaxt='n',xlab=expression(y[1]),ylab=expression(y[2]),cex.lab=1.25)
polygon(c(x,rev(x)),c(y2,rev(y1)),col="grey")
axis(side=1,at=c(0,2),labels=c("0",""))
axis(side=2,at=c(0,2),labels=c("0",""))
# Right
x <- seq(0,2,0.001)
y1 <- x
y2 <- -x
plot(x,y1,type="l",bty="n",xlim=c(-2,2),ylim=c(-2,2),xaxt='n',yaxt='n',xlab=expression(u[1]),ylab=expression(u[2]),cex.lab=1.25)
polygon(c(x,rev(x)),c(y2,rev(y1)),col="grey")
axis(side=1,at=c(-2,0,2),labels=c("","0",""))
axis(side=2,at=c(-2,0,2),labels=c("","0",""))


# Figure 6.14
# Page 38
y = seq(0,80,0.1)
pdf = dexp(y,1/12)
plot(y,pdf,type="l",xlab="y",ylab="PDF",cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)


# Figure 6.15
# Page 39
# Left
y = seq(0,6,0.1)
pdf = dexp(y,1/(12/14))
plot(y,pdf,type="l",xlab="y",ylab="PDF (minimum)",cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)
# Right
y = seq(0,120,0.1)
pdf = (14/12)*exp(-y/12)*(1-exp(-y/12))^13
plot(y,pdf,type="l",xlab="y",ylab="PDF (maximum)",cex.lab=1.25)
abline(h=0)


# Figure 6.16
# Page 41
# Upper left
y = seq(0,1,0.01)
pdf = rep(1,length(y))
plot(y,pdf,type="l",xlab="y",ylab="PDF",yaxt='n',xlim=c(0,1),ylim=c(0,3/2),cex.lab=1.25)
axis(side=2,at=c(0,1),labels=c("0","1"))
abline(h=0)
lines(c(1,1.5),c(0,0),lty=1)
# Other figures
n = 25 # number of observations
k = 13 # order statistic position
y = seq(0,1,0.001)
plot(y,dbeta(y,k,n-k+1),type="l",lty=1,xlab="y",ylab="PDF (median)",ylim=c(0,max(dbeta(y,k,n-k+1))),cex.lab=1.25)
abline(h=0)
# abline(v=0,lty=2) # use if k=1 (min)
# abline(v=1,lty=2) # use if k=n (max)


# Figure 6.17
# Page 43
# Left
x <- seq(0,1,0.001)
y1 <- x
y2 <- rep(1,length(x))
plot(x,y1,type="l",bty="n",xlim=c(0,1.1),ylim=c(0,1.1),xlab=expression(y[1]),ylab=expression(y[n]),xaxt='n',yaxt='n',cex.lab=1.25)
polygon(c(x,rev(x)),c(y2,rev(y1)),col="grey")
axis(side=1,at=c(0,1,1.2),labels=c("0","1",""))
axis(side=2,at=c(0,1,1.2),labels=c("0","1",""))
# Right
x <- seq(0,1,0.001)
y1 <- x
y2 <- rep(1,length(x))
plot(x,y1,type="l",bty="n",xlim=c(0,1.1),ylim=c(0,1.1),xlab=expression(u[1]),ylab=expression(u[2]),xaxt='n',yaxt='n',cex.lab=1.25)
polygon(c(x,rev(x)),c(y2,rev(y1)),col="grey")
axis(side=1,at=c(0,1,1.2),labels=c("0","1",""))
axis(side=2,at=c(0,1,1.2),labels=c("0","1",""))