###############################################
## Author: Joshua M. Tebbs
## Date: 11 August 2025
## Update: 17 August 2025
## STAT 509 course notes: R Code Chapter 6
###############################################

# 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 car R package
# See https://cran.r-project.org/web/packages/car/car.pdf for documentation

# Figure 6.1
# Page 89
# Left (PDF)
x = seq(40,280,0.01)
pdf = dnorm(x,160,30)
plot(x,pdf,type="l",xlab="x",ylab=expression(f[X](x)),xaxt="n",yaxt="n",cex.lab=1.25)
axis(1, at = c(160),labels = expression(mu),cex.axis=1.5)
abline(h=0)
lines(c(160,190),c(dnorm(190,160,30),dnorm(190,160,30)),lty=1)
text(175,dnorm(190,160,30)+0.00025,expression(sigma),cex=1.5)
# Right (Histogram)
strength = c(105,221,183,186,121,181,180,143,97,154,153,174,120,168,167,141,
	245,228,174,199,181,158,176,110,163,131,154,115,160,208,158,133,
	207,180,190,193,194,133,156,123,134,178,76,167,184,135,229,146,
	218,157,101,171,165,172,158,169,199,151,142,163,145,171,148,158,
	160,175,149,87,160,237,150,135,196,201,200,176,150,170,118,149)
bins = seq(0,300,20)
hist(strength,breaks=bins,xlab="Strength (in psi)",ylab="Count",ylim=c(0,25),main="",col="lightblue")


# Figure 6.2
# Page 92
# Left
y = seq(10,30,0.01)
plot(y,dnorm(y,20,2.5),type="l",lty=1,bty="n",xaxt="n",yaxt="n",xlab="",ylab="",ylim=c(-0.01,0.16),cex.lab=1.25)
abline(h=0)
points(x=20,y=0,pch=19,cex=1.25)
text(20,-0.01,expression(theta),cex=1.5)
text(30,-0.0075,expression(hat(theta)),cex=1.5)
# Right
y = seq(10,30,0.01)
plot(y,dnorm(y,20,2.5),type="l",lty=1,bty="n",xaxt="n",yaxt="n",xlab="",ylab="",ylim=c(-0.01,0.16),cex.lab=1.25)
abline(h=0)
points(x=17,y=0,pch=19,cex=1.25)
text(17,-0.01,expression(theta),cex=1.5)
text(30,-0.0075,expression(hat(theta)),cex=1.5)


# Figure 6.3
# Page 93
# Left
y = seq(10,30,0.01)
plot(y,dnorm(y,20,2.5),type="l",lty=1,bty="n",xaxt="n",yaxt="n",xlab="",ylab="",ylim=c(-0.01,0.16),cex.lab=1.25)
abline(h=0)
points(x=20,y=0,pch=19,cex=1.25)
text(20,-0.01,expression(theta),cex=1.5)
text(30,-0.0075,expression(hat(theta)),cex=1.5)
# Right
y = seq(10,30,0.01)
plot(y,dnorm(y,20,0.5),type="l",lty=1,bty="n",xaxt="n",yaxt="n",xlab="",ylab="",ylim=c(-0.04,0.8),cex.lab=1.25)
abline(h=0)
points(x=20,y=0,pch=19,cex=1.25)
text(20,-0.04,expression(theta),cex=1.5)
text(30,-0.035,expression(hat(theta)),cex=1.5)


# Figure 6.4
# Page 96
x = seq(24,88,0.1)
# Plot population distribution
plot(x,dnorm(x,56,8),type="l",lty=1,xlab="Ejection fraction",ylab="",ylim=c(0,0.2))
# Add sampling distribution
lines(x,dnorm(x,56,2),lty=4)
abline(h=0)
points(x=56,y=0,pch=19,cex=1.5)
# Add legend
legend(22.5,0.2,lty = c(1,4), 
	c(
	expression(paste("Population distribution")),
	expression(paste("Sampling distribution"))
	))


# Figure 6.5
# Page 98
# Left (Poisson PMF)
x = seq(0,10,1)
prob = dpois(x,1.5)
plot(x,prob,type="h",xlab="x",ylab=expression(p[X](x)),ylim=c(0,max(prob)),cex.lab=1.25,lwd=2)
abline(h=0)
points(x=1.5,y=0,pch=19,cex=1.5)
# Right (Sampling distribution)
x = seq(0.8,2.2,0.0001)
pdf = dnorm(x,1.5,sqrt(0.03))
plot(x,pdf,type="l",xlab="Average number of defects (per square meter)",ylab="",cex.lab=1.25)
abline(h=0)
x = seq(2,2.2,0.0001)
y = dnorm(x,1.5,sqrt(0.03))
polygon(c(2,x,2.2),c(0,y,0),col="lightblue")
points(x=2,y=0,pch=19,cex=1)
text(2.1,0.1,0.002,cex=1)


# Figure 6.6
# Page 99
x = seq(0,120,0.1)
pdf = dexp(x,0.05)
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)
points(x=20,y=0,pch=19,cex=1.5)


# Figure 6.7
# Page 100
x = seq(20-4*sqrt(10),20+4*sqrt(10),0.01)
pdf = dnorm(x,20,sqrt(10))
plot(x,pdf,type="l",xlab="Average time to enrollment (in days)",ylab="",cex.lab=1.25)
abline(h=0)


# Figure 6.8
# Page 102
t = seq(-5,5,0.01)
plot(t,dnorm(t,0,1),type="l",lty=1,xlab="t",ylab=expression(f[T](t)),cex.lab=1.25)
lines(t,dt(t,3),lty=4)
lines(t,dt(t,10),lty=8)
abline(h=0)
# Add legend
legend(2,0.35,lty=c(1,4,8),
	c(expression(paste("N(0,1)")),
	expression(paste("t(3)")),
	expression(paste("t(10)"))
	))


# Figure 6.9
# Page 104
# QQ plot for strength data in Example 6.1
strength = c(105,221,183,186,121,181,180,143,97,154,153,174,120,168,167,141,
	245,228,174,199,181,158,176,110,163,131,154,115,160,208,158,133,
	207,180,190,193,194,133,156,123,134,178,76,167,184,135,229,146,
	218,157,101,171,165,172,158,169,199,151,142,163,145,171,148,158,
	160,175,149,87,160,237,150,135,196,201,200,176,150,170,118,149)
qqPlot(strength,distribution="norm",mean=mean(strength),sd=sd(strength),
	xlab="Normal quantiles",ylab="Strength (in psi)",pch=16,
	envelope=list(border=TRUE,style="lines"),id=FALSE)


# Figure 6.10
# Page 105
arsenic = c(17.6,10.4,13.5,4,19.9,16,12,12.2,11.4,12.7,3,10.3,21.4,19.4,9,6.5,10.1,8.7,9.7,6.4,
	9.7,63,15.5,10.7,18.2,7.5,6.1,6.7,6.9,0.8,73.5,12,28,12.6,9.4,6.2,15.3,7.3,10.7,15.9,
	5.8,1,8.6,1.3,13.7,2.8,2.4,1.4,2.9,13.1,15.3,9.2,11.7,4.5,1,1.2,0.8,1,2.4,4.4,2.2,2.9,
	3.6,2.5,1.8,5.9,2.8,1.7,4.6,5.4,3,3.1,1.3,2.6,1.4,2.3,1,5.4,1.8,2.6,3.4,1.4,10.7,18.2,
	7.7,6.5,12.2,10.1,6.4,10.7,6.1,0.8,12,28.1,9.4,6.2,7.3,9.7,62.1,15.5,6.4,9.5)
# Left (Histogram)
bins = seq(0,80,5)
hist(arsenic,breaks = bins,xlab="Arsenic concentration (in ppb)",ylab="Count",main="",col="lightblue")
# Right (QQ plot)
qqPlot(arsenic,distribution="norm",mean=mean(arsenic),sd=sd(arsenic),
	xlab="Normal quantiles",ylab="Arsenic concentration (in ppb)",pch=16,
	envelope=list(border=TRUE,style="lines"),id=FALSE)