###############################################
## Author: Joshua M. Tebbs
## Date: 11 August 2026
## Update: 29 August 2026
## STAT 515 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 "car" R package
# For normal qq plots
# See https://cran.r-project.org/web/packages/car/car.pdf for documentation

# Figure 5.1
# Page 53
x = seq(-5,5,0.001)
pdf = dt(x,10)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(-5,qt(0.025,10),0.001)
y = dt(x,10)
polygon(c(-5,x,qt(0.025,10)),c(0,y,0),col="lightblue")
points(x=qt(0.025,10),y=0,pch=19,cex=1)
x = seq(qt(0.975,10),5,0.001)
y = dt(x,10)
polygon(c(qt(0.975,10),x,5),c(0,y,0),col="lightblue")
points(x=qt(0.975,10),y=0,pch=19,cex=1)
text(-0.025,0.1,expression(1-alpha),cex=1.25)
text(-3.5,0.05,expression(alpha/2),cex=1.25)
text(3.5,0.05,expression(alpha/2),cex=1.25)
text(1.5,0.3,"N(0,1)",cex=1.25)


# Example 5.1
# Page 54-55
bricks = c(4.54,4.64,4.58,4.78,4.58,4.62,4.55,4.63,4.51,4.49,4.50,4.51,4.63,4.47,4.36,
	4.61,4.53,4.45,4.26,4.40,4.48,4.63,4.47,4.46,4.57,4.41,4.50,4.62,4.50,4.61,4.49,
	4.79,4.39,4.70,4.39,4.45)
options(digits=3)
mean(bricks) # sample mean
sd(bricks) # sample standard deviation

# N(0,1) pdf with 95% confidence
# Page 55
x = seq(-3.5,3.5,0.001)
pdf = dnorm(x,0,1)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(-3.5,qnorm(0.025,0,1),0.001)
y = dnorm(x,0,1)
polygon(c(-3.5,x,qnorm(0.025,0,1)),c(0,y,0),col="lightblue")
points(x=qnorm(0.025,0,1),y=0,pch=19,cex=1)
x = seq(qnorm(0.975,0,1),3.5,0.001)
y = dnorm(x,0,1)
polygon(c(qnorm(0.975,0,1),x,3.5),c(0,y,0),col="lightblue")
points(x=qnorm(0.975,0,1),y=0,pch=19,cex=1)
text(-0.025,0.1,0.95,cex=1.25)
text(-3,0.04,0.025,cex=1.25)
text(3,0.04,0.025,cex=1.25)
text(1.5,0.3,"N(0,1)",cex=1.25)
text(qnorm(0.975,0,1),-0.011,1.96,cex=1.1)
text(qnorm(0.025,0,1),-0.011,-1.96,cex=1.1)


# Figure 5.2
# Page 58
bricks = c(4.54,4.64,4.58,4.78,4.58,4.62,4.55,4.63,4.51,4.49,4.50,4.51,4.63,4.47,4.36,
	4.61,4.53,4.45,4.26,4.40,4.48,4.63,4.47,4.46,4.57,4.41,4.50,4.62,4.50,4.61,4.49,
	4.79,4.39,4.70,4.39,4.45)
# Left (Histogram)
bins = seq(4,5,0.1)
hist(bricks,breaks=bins,xlab="Weight (in lbs)",ylab="Count",ylim=c(0,12),xlim=c(4,5),main="",col="lightblue")
# Right (QQ plot)
qqPlot(bricks,distribution="norm",mean=mean(bricks),sd=sd(bricks),
	xlab="Normal quantiles",ylab="Weight (in lbs)",pch=16,
	envelope=list(border=TRUE,style="lines"),id=FALSE)


# Figure 5.3
# Page 62
# This is the same figure as Figure 5.1.


# Example 5.3
# Page 63
cholesterol = c(211,228,284,95,181,122,202,215,224,136,235,133,168,181,169,142,181,120,240,102,
		    135,227,211,197,175,209,186,166,208,129,181,226,204,111,267,131,140,163,187,
		    175,208,176,85,225,211,192,152,209,172,132,147,224,191,232,100,163,130,224,156,
		    294,181,171,149,206,271,187,151,82,206,176,166,168,228,178,156,124,87,127,168,
		    243,165,209,176,131,136,171,187,191,125,222,212,141,174,95,154,133)
mean(cholesterol) # sample mean
# N(0,1) pdf with alpha = 0.05 rejection region
# The same figure is on page 55.


# Figure 5.4
# Page 65
# Left
x = seq(-5,5,0.001)
pdf = dt(x,10)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(-5,qt(0.05,10),0.001)
y = dt(x,10)
polygon(c(-5,x,qt(0.05,10)),c(0,y,0),col="lightblue")
points(x=qt(0.05,10),y=0,pch=19,cex=1)
text(-0.025,0.1,expression(1-alpha),cex=1.25)
text(-3,0.05,expression(alpha),cex=1.25)
text(1.5,0.3,"N(0,1)",cex=1.25)
# Right
x = seq(-5,5,0.001)
pdf = dt(x,10)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(qt(0.95,10),5,0.001)
y = dt(x,10)
polygon(c(qt(0.95,10),x,5),c(0,y,0),col="lightblue")
points(x=qt(0.95,10),y=0,pch=19,cex=1)
text(-0.025,0.1,expression(1-alpha),cex=1.25)
text(3,0.05,expression(alpha),cex=1.25)
text(1.5,0.3,"N(0,1)",cex=1.25)


# Example 5.4
# Page 65-66
times = c(2.7,2.4,3.4,3.2,2.8,3.9,3.6,2.8,2.9,2.2,3.3,3.6,4.5,3.1,3.6,3.4,3.1,3.0,2.1,2.5,3.8,3.3,3.2,3.1,2.7,
	    2.1,2.9,3.5,2.9,2.7,3.8,3.1,2.8,1.2,3.6,2.3,2.7,2.5,2.3,2.2,2.2,3.0,2.5,2.4,2.1,2.7,2.9,3.6,2.7,3.1)
mean (times) # sample mean

# One-sided rejection region (lower tail)
# Page 66
x = seq(-3.5,3.5,0.001)
pdf = dnorm(x,0,1)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(-3.5,qnorm(0.1,0,1),0.001)
y = dnorm(x,0,1)
polygon(c(-3.5,x,qnorm(0.1,0,1)),c(0,y,0),col="lightblue")
points(x=qnorm(0.1,0,1),y=0,pch=19,cex=1)
text(0.025,0.1,"0.90",cex=1.25)
text(-2.5,0.06,"0.10",cex=1.25)
text(1.5,0.3,"N(0,1)",cex=1.25)
text(qnorm(0.1,0,1),-0.011,-1.28,cex=1.1)


# Example 5.5
# Page 67-68
indent = c(22.4,23.6,24.0,24.9,25.5,25.6,25.8,26.1,26.4,26.7,27.4,27.6,28.3,
	   29.0,29.1,29.6,29.7,29.8,29.9,30.0,30.4,30.5,30.7,30.7,31.0,31.0,
	   31.4,31.6,31.7,31.9,31.9,32.0,32.1,32.4,32.5,32.5,32.6,32.9,33.1,
	   33.3,33.5,33.5,33.5,33.5,33.6,33.6,33.8,33.9,34.1,34.2,34.6,34.6,
	   35.0,35.2,35.2,35.4,35.4,35.4,35.5,35.7,35.8,36.0,36.0,36.0,36.1,
	   36.1,36.2,36.4,36.6,37.0,37.4,37.5,37.5,38.0,38.7,38.8,39.8,41.0,
	   42.0,42.1,44.6,48.3,55.0)
mean(indent) # sample mean
sd(indent) # sample standard deviation

# One-sided rejection region (upper tail)
# Page 68
x = seq(-3.5,3.5,0.001)
pdf = dnorm(x,0,1)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(qnorm(0.99,0,1),3.5,0.001)
y = dnorm(x,0,1)
polygon(c(qnorm(0.99,0,1),x,3.5),c(0,y,0),col="lightblue")
points(x=qnorm(0.99,0,1),y=0,pch=19,cex=1)
text(0.025,0.1,"0.99",cex=1.25)
text(3.1,0.03,"0.01",cex=1.25)
text(1.5,0.3,"N(0,1)",cex=1.25)
text(qnorm(0.99,0,1),-0.011,2.33,cex=1.1)


# Figure 5.5
# Page 69
# Left (Histogram)
bins = seq(20,60,2)
hist(indent,breaks=bins,xlab="Indentation (in mm)",ylab="Count",xlim=c(20,60),main="",col="lightblue")
# Right (QQ plot)
library(car)
qqPlot(indent,distribution="norm",mean=mean(indent),sd=sd(indent),
    xlab="Normal quantiles",ylab="Indentation (in mm)",pch=16,
    envelope=list(border=TRUE,style="lines"),id=FALSE)


# Figure 5.6
# Page 71
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 5.7
# Page 72
x = seq(-5,5,0.001)
pdf = dt(x,10)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(-5,qt(0.025,10),0.001)
y = dt(x,10)
polygon(c(-5,x,qt(0.025,10)),c(0,y,0),col="lightblue")
points(x=qt(0.025,10),y=0,pch=19,cex=1)
x = seq(qt(0.975,10),5,0.001)
y = dt(x,10)
polygon(c(qt(0.975,10),x,5),c(0,y,0),col="lightblue")
points(x=qt(0.975,10),y=0,pch=19,cex=1)
text(-0.025,0.1,expression(1-alpha),cex=1.25)
text(-3.5,0.05,expression(alpha/2),cex=1.25)
text(3.5,0.05,expression(alpha/2),cex=1.25)
text(1.5,0.3,"t(n-1)",cex=1.25)


# Example 5.6
# Page 73
twins = c(20.7,16.7,22.5,12.1,2.9,23.5,6.4,1.3,39.6,35.6,
	  10.9,7.1,46.0,23.4,29.4,44.1,13.8,24.3,9.3,3.4)
mean(twins) # sample mean
sd(twins) # sample standard deviation


# Figure 5.8
# Page 75
library(car)
qqPlot(twins,distribution="norm",mean=mean(twins),sd=sd(twins),
    xlab="Normal quantiles",ylab="Attention times (in hours)",pch=16,
    envelope=list(border=TRUE,style="lines"),id=FALSE)


# Figure 5.9
# Page 76
# Top--this is the same as Figure 5.7.
# Bottom left
x = seq(-5,5,0.001)
pdf = dt(x,10)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(-5,qt(0.05,10),0.001)
y = dt(x,10)
polygon(c(-5,x,qt(0.05,10)),c(0,y,0),col="lightblue")
points(x=qt(0.05,10),y=0,pch=19,cex=1)
text(-0.025,0.1,expression(1-alpha),cex=1.25)
text(-3,0.05,expression(alpha),cex=1.25)
text(1.5,0.3,"t(n-1)",cex=1.25)
# Bottom right
x = seq(-5,5,0.001)
pdf = dt(x,10)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(qt(0.95,10),5,0.001)
y = dt(x,10)
polygon(c(qt(0.95,10),x,5),c(0,y,0),col="lightblue")
points(x=qt(0.95,10),y=0,pch=19,cex=1)
text(-0.025,0.1,expression(1-alpha),cex=1.25)
text(3,0.05,expression(alpha),cex=1.25)
text(1.5,0.3,"t(n-1)",cex=1.25)


# Example 5.7
# Page 76-77
cycles = c(2226,2283,875,733,1390,1744,1174,1468,1229,1386,
	     1843,914,1518,1288,1507,1321,819,1370,1543,2986)
mean(cycles) # sample mean
sd(cycles) # sample standard deviation

# One-sided rejection region (upper tail) using alpha = 0.05
# Page 77
x = seq(-4,4,0.001)
pdf = dt(x,19)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(qt(0.95,19),4,0.001)
y = dt(x,19)
polygon(c(qt(0.95,19),x,4),c(0,y,0),col="lightblue")
points(x=qt(0.95,19),y=0,pch=19,cex=1)
text(-0.025,0.1,"0.95",cex=1.25)
text(3,0.05,"0.05",cex=1.25)
text(1.5,0.3,"t(19)",cex=1.25)
text(qt(0.95,19),-0.011,1.729,cex=1.1)


# Figure 5.10
# Page 78
library(car)
qqPlot(cycles,distribution="norm",mean=mean(cycles),sd=sd(cycles),
    xlab="Normal quantiles",ylab="Number of cycles",pch=16,
    envelope=list(border=TRUE,style="lines"),id=FALSE)

# One-sided rejection region (upper tail) using alpha = 0.01
# Page 79
x = seq(-4,4,0.001)
pdf = dt(x,19)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(qt(0.99,19),4,0.001)
y = dt(x,19)
polygon(c(qt(0.99,19),x,4),c(0,y,0),col="lightblue")
points(x=qt(0.99,19),y=0,pch=19,cex=1)
text(0.025,0.1,"0.99",cex=1.25)
text(3.2,0.03,"0.01",cex=1.25)
text(1.5,0.3,"t(19)",cex=1.25)
text(qt(0.99,19),-0.011,2.539,cex=1.1)


# Figure 5.11
# Page 82
x = seq(-3.5,3.5,0.001)
pdf = dnorm(x,0,1)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",xlim = c(-3.5,6),ylim=c(0,0.4))
abline(h=0)
x = seq(qnorm(0.95,0,1),3.5,0.001)
y = dnorm(x,0,1)
polygon(c(qnorm(0.95,0,1),x,3.5),c(0,y,0),col="lightblue")
points(x=qnorm(0.95,0,1),y=0,pch=19,cex=1)
text(2.5,0.035,expression(alpha),cex=1.25)
text(-1.25,0.3,expression(H[0]),cex=1.25)
x = seq(-2.5,5.5,0.001)
lines(x,dnorm(x,2,1),lty=4)
text(3.25,0.3,expression(H[a]),cex=1.25)
axis(side=1,at=c(0,2),labels=c(expression(mu[0]),expression(mu[a])),cex.axis=1.25)


# Figure 5.12
# Page 84
mu = seq(12,22,0.001)
mu.0 = 15
n = 25
alpha = 0.05
sigma = 5
plot(mu,1-pnorm(qnorm(1-alpha) + (mu.0-mu)/(sigma/sqrt(n))),
	type="l",lty=1,xlab=expression(mu),ylab=expression(PWR(mu)),xaxp=c(12,22,10),ylim=c(0,1),cex.lab=1.25)
abline(h=alpha,lty=2)
points(x=18,y=1-pnorm(qnorm(1-alpha) + (mu.0-18)/(sigma/sqrt(n))),pch=19,cex=1)


# Example 5.10
# Page 87-88
salmonella = c(0.593,0.142,0.329,0.691,0.231,0.793,0.519,0.392,0.418)
mean(salmonella) # sample mean
sd(salmonella) # sample standard deviation
t.test(salmonella, alternative="greater", mu = 0.3, conf.level = 0.95)

# Page 88
# Left: One-sided rejection region (upper tail) using alpha = 0.05
x = seq(-4,4,0.001)
pdf = dt(x,8)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(qt(0.95,8),4,0.001)
y = dt(x,8)
polygon(c(qt(0.95,8),x,4),c(0,y,0),col="lightblue")
points(x=qt(0.95,8),y=0,pch=19,cex=1)
text(-0.025,0.1,"0.95",cex=1.25)
text(3,0.05,"0.05",cex=1.25)
text(1.5,0.3,"t(8)",cex=1.25)
text(qt(0.95,8),-0.011,"1.860",cex=1.1)
# Right: One-sided rejection region (upper tail) using alpha = 0.01
x = seq(-4,4,0.001)
pdf = dt(x,8)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(qt(0.99,8),4,0.001)
y = dt(x,8)
polygon(c(qt(0.99,8),x,4),c(0,y,0),col="lightblue")
points(x=qt(0.99,8),y=0,pch=19,cex=1)
text(-0.025,0.1,"0.99",cex=1.25)
text(3.5,0.025,"0.01",cex=1.25)
text(1.5,0.3,"t(8)",cex=1.25)
text(qt(0.99,8),-0.011,2.897,cex=1.1)


# Figure 5.13
# Page 89
x = seq(-4,4,0.001)
pdf = dt(x,8)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(2.2,4,0.001)
y = dt(x,8)
polygon(c(2.2,x,4),c(0,y,0),col="lightblue")
points(x=2.2,y=0,pch=19,cex=1)
text(3.3,0.035,"0.0295",cex=1.25)
text(1.5,0.3,"t(8)",cex=1.25)
text(2.2,-0.011,"2.20",cex=1.1)