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

# 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 

# Figure 4.1
# Page 49
x = seq(0,30,0.01)
pdf = 3000/(10+x)^4
plot(x,pdf,type="l",xlab="x",ylab=expression(f[X](x)),cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)
# Add shaded area
x = seq(10,15,0.001)
y = 3000/(10+x)^4
polygon(c(10,x,15),c(0,y,0),col="lightblue")
points(x=c(10,15),y=c(0,0),pch=19,cex=1)


# Figure 4.2
# Page 50
x = seq(4.9,5.1,0.01)
pdf = 5+x-x
plot(x,pdf,type="l",xlab="x",ylab=expression(f[X](x)),xaxp=c(4.9,5.1,4),xlim=c(4.875,5.125),yaxp=c(5,5,1),
	ylim=c(0,12.5),cex.lab=1.25)
abline(h=0,lty=2)
# Add shaded area
x = seq(4.9,4.95,0.001)
y = 5+x-x
polygon(c(4.9,x,4.95),c(0,y,0),col="lightblue")
points(x=c(4.9,4.95),y=c(0,0),pch=19,cex=1)


# Figure 4.3
# Page 52
x = seq(0,1,0.001)
pdf = 90*x^8*(1-x)
plot(x,pdf,type="l",xlab="x",ylab=expression(f[X](x)),xaxp=c(0,1,4),xlim=c(0,1),cex.lab=1.25)
abline(h=0)
# Add shaded area
x = seq(0.9,1,0.001)
y = 90*x^8*(1-x)
polygon(c(0.9,x,1),c(0,y,0),col="lightblue")
points(x=c(0.9,1),y=c(0,0),pch=19,cex=1)


# Figure 4.4
# Page 54
x = seq(0,30,0.01)
# Left (PDF)
pdf = 3000/(10+x)^4
plot(x,pdf,type="l",xlab="x",ylab=expression(f[X](x)),cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)
# Right (CDF)
cdf = 1-(10/(10+x))^3
plot(x,cdf,type="l",xlab="x",ylab=expression(F[X](x)),ylim=c(0,1),cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)


# Figure 4.5
# Page 55
x = seq(4.9,5.1,0.01)
# Left (PDF)
pdf = 5+x-x
plot(x,pdf,type="l",xlab="x",ylab=expression(f[X](x)),xaxp=c(4.9,5.1,4),xlim=c(4.875,5.125),yaxp=c(5,5,1),
	ylim=c(0,12.5),cex.lab=1.25)
abline(h=0,lty=2)
# Right (CDF)
cdf = 5*x-24.5
plot(x,cdf,type="l",xlab="x",ylab=expression(F[X](x)),xaxp=c(4.9,5.1,4),xlim=c(4.875,5.125),
	ylim=c(0,1.33),yaxp=c(0,1,4),cex.lab=1.25)
lines(c(5.1,5.2),c(1,1),lty=1)
abline(h=0)


# Figure 4.6
# Page 55
x = seq(0,1,0.001)
# Left (PDF)
pdf = 90*x^8*(1-x)
plot(x,pdf,type="l",xlab="x",ylab=expression(f[X](x)),xaxp=c(0,1,4),xlim=c(0,1),cex.lab=1.25)
abline(h=0)
# Right (CDF)
cdf = 10*x^9-9*x^(10)
plot(x,cdf,type="l",xlab="x",ylab=expression(F[X](x)),xaxp=c(0,1,4),xlim=c(0,1.15),ylim=c(0,1),cex.lab=1.25)
lines(c(1,1.15),c(1,1),lty=1)
abline(h=0)


# Figure 4.7
# Page 57
x = seq(12.5,13,0.001)
pdf = 20*exp(-20*(x-12.5))
plot(x,pdf,type="l",xlab="x",ylab=expression(f[X](x)),cex.lab=1.25)
abline(h=0)
abline(v=12.5,lty=2)
points(x=12.55,y=0,pch=19,cex=1.5)


# Figure 4.8
# Page 59
x = seq(400,600,0.01)
pdf = 1200/x^2
plot(x,pdf,type="l",xlab="x",ylab=expression(f[X](x)),xlim=c(390,610),ylim=c(0,0.01),cex.lab=1.25)
abline(h=0)
lines(c(400,400),c(0,1200/400^2),lty=2)
lines(c(600,600),c(0,1200/600^2),lty=2)
points(x=486.6,y=0,pch=19,cex=1.5)


# Figure 4.9
# Page 61
x = seq(12.5,13,0.001)
# Left (PDF)
pdf = 20*exp(-20*(x-12.5))
plot(x,pdf,type="l",xlab="x",ylab=expression(f[X](x)),cex.lab=1.25)
abline(h=0)
abline(v=12.5,lty=2)
x = seq(12.5,12.5-log(0.05)/20,0.001)
y = 20*exp(-20*(x-12.5))
polygon(c(12.5,x,12.5-log(0.05)/20),c(0,y,0),col="lightblue")
points(x=12.5-log(0.05)/20,y=0,pch=19,cex=1)
text(12.56,2,0.95,cex=1.25)
# Right (CDF)
x = seq(12.5,13,0.001)
cdf = 1-exp(-20*(x-12.5))
plot(x,cdf,type="l",xlab="x",ylab=expression(F[X](x)),xlim=c(12.45,12.9),ylim=c(0,1),cex.lab=1.25)
abline(h=0)
points(x=12.5-log(0.05)/20,y=0,pch=19,cex=1)
lines(c(12.5-log(0.05)/20,12.5-log(0.05)/20),c(0,0.95),lty=2)
lines(c(12.4,12.5-log(0.05)/20),c(0.95,0.95),lty=2)
text(12.45,0.935,0.95,cex=1)


# Figure 4.10
# Page 62
x = seq(0,8,0.01)
plot(x,dexp(x,2),type="l",lty=1,xlab="x",ylab=expression(f[X](x)),cex.lab=1.25)
lines(y,dexp(y,1),lty=4)
lines(y,dexp(y,1/2),lty=8)
abline(h=0)
abline(v=0,lty=2)
# Add legend
legend(3,1,lty = c(1,4,8), 
	c(
	expression(paste(lambda, "=2")),
	expression(paste(lambda, "=1")),
	expression(paste(lambda, "=0.5"))
	))


# Figure 4.11
# Page 63
x = seq(0,25,0.01)
# Left (PDF)
pdf = dexp(x,0.25)
plot(x,pdf,type="l",lty=1,xlab="x",ylab=expression(f[X](x)),cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)
x = seq(10,25,0.01)
y = dexp(x,0.25)
polygon(c(10,x,25),c(0,y,0),col="lightblue")
points(x=10,y=0,pch=19,cex=1)
text(15,0.02,0.082,cex=1.25)
# Right (CDF)
x = seq(0,25,0.01)
cdf = pexp(x,0.25)
plot(x,cdf,type="l",xlab="x",ylab=expression(F[X](x)),xlim=c(0,25),ylim=c(0,1),cex.lab=1.25)
abline(h=0)
points(x=2.77,y=0,pch=19,cex=1)
lines(c(2.77,2.77),c(0,0.50),lty=2)
lines(c(-1,2.77),c(0.50,0.50),lty=2)
text(-0.2,0.52,0.50,cex=1)


# Figure 4.12
# Page 65
x = seq(0,0.4,0.001)
# Left (PDF)
pdf = dexp(x,25)
plot(x,pdf,type="l",xlab="x",ylab=expression(f[X](x)),cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)
x = seq(1/6,0.4,0.001)
y = dexp(x,25)
polygon(c(1/6,x,0.5),c(0,y,0),col="lightblue")
points(x=1/6,y=0,pch=19,cex=1)
text(0.225,1.5,0.016,cex=1)


# Figure 4.13
# Page 67
x = seq(0,25,0.01)
plot(x,dgamma(x,1.5,0.6),type="l",lty=1,xlab="x",ylab=expression(f[X](x)),cex.lab=1.25)
lines(x,dgamma(x,2,0.5),lty=4)
lines(x,dgamma(x,2.5,0.4),lty=8)
abline(h=0)
# Add legend
legend(10,0.20,lty = c(1,4,8), 
	c(
	expression(paste(r, "=1.5, ", lambda, "=0.6")),
	expression(paste(r, "=2.0, ", lambda, "=0.5")),
	expression(paste(r, "=2.5, ", lambda, "=0.4"))
	))


# Figure 4.14
# Page 68
x = seq(0,25,0.01)
plot(x,dgamma(x,2.2,0.4),type="l",lty=1,xlab="x",ylab=expression(f[X](x)),cex.lab=1.25)
abline(h=0)


# Figure 4.15
# Page 70
x = seq(-10,10,0.01)
plot(x,dnorm(x,0,1),type="l",lty=1,xlab="x",ylab=expression(f[X](x)),cex.lab=1.25)
lines(x,dnorm(x,-2,2),lty=4)
lines(x,dnorm(x,1,3),lty=8)
abline(h=0)
# Add legend
legend(3.5,0.30,lty = c(1,4,8), 
	c(
	expression(paste(mu, "=0, ", sigma, "=1")),
	expression(paste(mu, "=-2, ", sigma, "=2")),
	expression(paste(mu, "=1, ", sigma, "=3"))
	))


# Figure 4.16
# Page 72
x = seq(319,575,0.05)
pdf = dnorm(x,447,32)
plot(x,pdf,type="l",xlab="x",ylab=expression(f[X](x)),cex.lab=1.25)
abline(h=0)
x = seq(319,400,0.05)
y = dnorm(x,447,32)
polygon(c(319,x,400),c(0,y,0),col="lightblue")
points(x=400,y=0,pch=19,cex=1)
text(360,0.0015,0.071,cex=1)