###############################################
## Author: Joshua M. Tebbs
## Date: 3 October 2026
## Update: 6 October 2026
## STAT 515 course notes: R Code Chapter 7
###############################################

# 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 7.1
# Page 125
q = seq(0,45,0.01)
plot(q,dchisq(q,5),type="l",lty=1,xlab=expression(chi^2),ylab=expression(f(chi^2)),ylim=c(0,0.16))
lines(q,dchisq(q,10),lty=4)
lines(q,dchisq(q,20),lty=8)
abline(h=0)
# Add legend
legend(25,0.125,lty = c(1,4,8), 
	c(
	expression(paste(nu, " = 5")),
	expression(paste(nu, " = 10")),
	expression(paste(nu, " = 20"))
	))


# Figure 7.2
# Page 126
x = seq(0,30,0.01)
pdf = dchisq(x,10)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",xaxt="n",yaxt="n",bty="n",ylim=c(0,0.1))
abline(h=0)
x = seq(0,qchisq(0.025,10),0.01)
y = dchisq(x,10)
polygon(c(0,x,qchisq(0.025,10)),c(0,y,0),col="lightblue")
points(x=qchisq(0.025,10),y=0,pch=19,cex=1)
x = seq(qchisq(0.975,10),30,0.01)
y = dchisq(x,10)
polygon(c(qchisq(0.975,10),x,30),c(0,y,0),col="lightblue")
points(x=qchisq(0.975,10),y=0,pch=19,cex=1)
text(9.8,0.02,expression(1-alpha),cex=1.25)
text(0.25,0.01,expression(alpha/2),cex=1.25)
text(25,0.0075,expression(alpha/2),cex=1.25)
text(15,0.08,expression(paste(chi^2, "(n-1)")),cex=1.25)


# Example 7.1
# Page 127-130
glucose = c(154.5,168.1,197.1,168.2,162.1,162.6,155.5,164.6,184.4,174.1, 
		160.9,171.3,213.4,155.2,180.2,196.6,197.0,164.7,171.4,208.6)
var(glucose) # sample variance

# chi^2(19) pdf with 95% confidence
# Page 128
x = seq(0,45,0.001)
pdf = dchisq(x,19)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",xaxt="n",yaxt="n",bty="n")
abline(h=0)
x = seq(0,qchisq(0.025,19),0.001)
y = dchisq(x,19)
polygon(c(0,x,qchisq(0.025,19)),c(0,y,0),col="lightblue")
points(x=qchisq(0.025,19),y=0,pch=19,cex=1)
x = seq(qchisq(0.975,19),45,0.001)
y = dchisq(x,19)
polygon(c(qchisq(0.975,19),x,45),c(0,y,0),col="lightblue")
points(x=qchisq(0.975,19),y=0,pch=19,cex=1)
text(26,0.05,expression(paste(chi^2, "(19)")),cex=1.25)
text(19.5,0.0175,0.95,cex=1.25)
text(3,0.008,0.025,cex=1.25)
text(40,0.006,0.025,cex=1.25)
text(qchisq(0.025,19),-0.0015,8.907,cex=1)
text(qchisq(0.975,19),-0.0015,32.85,cex=1)

# Figure 7.3
# Page 129
library(car)
qqPlot(glucose,distribution="norm",mean=mean(glucose),sd=sd(glucose),
	xlab="Normal quantiles",ylab="Blood glucose level (in mg/dL)",pch=16,
	envelope=list(border=TRUE,style="lines"),id=FALSE)

# Population variance confidence interval function I wrote:
# Page 130
var.ci = function(data,conf.level=0.95){
    df = length(data)-1
    chi.lower = qchisq((1-conf.level)/2,df)
    chi.upper = qchisq((1+conf.level)/2,df)
    s2 = var(data)
    c(df*s2/chi.upper,df*s2/chi.lower)
    }
var.ci(glucose,conf.level=0.95) # CI for population variance
sqrt(var.ci(glucose,conf.level=0.95)) # CI for population standard deviation


# Figure 7.4
# Page 131
# Top--this is the same as Figure 7.2.
# Bottom left
x = seq(0,30,0.01)
pdf = dchisq(x,10)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",xaxt="n",yaxt="n",bty="n",ylim=c(0,0.1))
abline(h=0)
x = seq(0,qchisq(0.05,10),0.01)
y = dchisq(x,10)
polygon(c(0,x,qchisq(0.05,10)),c(0,y,0),col="lightblue")
points(x=qchisq(0.05,10),y=0,pch=19,cex=1)
text(9.8,0.02,expression(1-alpha),cex=1.25)
text(0.25,0.01,expression(alpha),cex=1.25)
text(15,0.08,expression(paste(chi^2, "(n-1)")),cex=1.25)
# Bottom right
x = seq(0,30,0.01)
pdf = dchisq(x,10)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",xaxt="n",yaxt="n",bty="n",ylim=c(0,0.1))
abline(h=0)
x = seq(qchisq(0.95,10),30,0.01)
y = dchisq(x,10)
polygon(c(qchisq(0.95,10),x,30),c(0,y,0),col="lightblue")
points(x=qchisq(0.95,10),y=0,pch=19,cex=1)
text(9.8,0.02,expression(1-alpha),cex=1.25)
text(25,0.0075,expression(alpha),cex=1.25)
text(15,0.08,expression(paste(chi^2, "(n-1)")),cex=1.25)


# Example 7.2
# Page 131-134
concentration = c(91.28,92.83,89.35,91.90,82.85,94.83,89.83,89.00,84.62,86.96, 
			88.32,91.17,83.86,89.74,92.24,92.59,84.21,89.36,90.96,92.85,
			89.39,89.82,89.91,92.16,88.67,89.35,86.51,89.04,91.82,93.02,
			88.32,88.76,89.26,90.36,87.16,91.74,86.12,92.10,83.33,87.61,
			88.20,92.78,86.35,93.84,91.20,93.44,86.77,83.77,93.19,81.79)
var(concentration) # sample variance

# One-sided rejection region (lower tail) using alpha = 0.05
# Page 132
x = seq(0,100,0.001)
pdf = dchisq(x,49)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",xaxt="n",yaxt="n",bty="n")
abline(h=0)
x = seq(0,qchisq(0.05,49),0.001)
y = dchisq(x,49)
polygon(c(0,x,qchisq(0.05,49)),c(0,y,0),col="lightblue")
points(x=qchisq(0.05,49),y=0,pch=19,cex=1)
text(66,0.03,expression(paste(chi^2, "(49)")),cex=1.25)
text(50,0.009,0.95,cex=1.25)
text(20,0.004,0.05,cex=1.25)
text(qchisq(0.05,49),-0.001,33.93,cex=1)

# P-value
# Page 133
x = seq(0,100,0.001)
pdf = dchisq(x,49)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",xaxt="n",yaxt="n",bty="n")
abline(h=0)
x = seq(0,55.18,0.001)
y = dchisq(x,49)
polygon(c(0,x,55.18),c(0,y,0),col="lightblue")
points(x=55.18,y=0,pch=19,cex=1)
text(66,0.03,expression(paste(chi^2, "(49)")),cex=1.25)
text(55.18,-0.001,55.18,cex=1)
text(43,0.01,0.747,cex=1.25)

# Figure 7.5
# Page 134
library(car)
qqPlot(concentration,distribution="norm",mean=mean(concentration),sd=sd(concentration),
	xlab="Normal quantiles",ylab="Drug concentration (measured as a %)",pch=16,
	envelope=list(border=TRUE,style="lines"),id=FALSE)


# Figure 7.6
# Page 135
q = seq(0,6,0.01)
plot(q,df(q,5,10),type="l",lty=1,xlab="f",ylab="F pdf",cex.lab=1.25,ylim=c(0,1.1))
lines(q,df(q,10,20),lty=4)
lines(q,df(q,20,40),lty=8)
abline(h=0)
# Add legend
legend(3,0.8,lty = c(1,4,8),
	c(expression(paste(nu[1]," = 5, ",nu[2]," = 10")),
	expression(paste(nu[1]," = 10, ",nu[2], " = 20")),
	expression(paste(nu[1]," = 20, ", nu[2]," = 40"))
	))


# Figure 7.7
# Page 136
x = seq(0,30,0.01)
pdf = dchisq(x,10)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",xaxt="n",yaxt="n",bty="n",ylim=c(0,0.1))
abline(h=0)
x = seq(0,qchisq(0.025,10),0.01)
y = dchisq(x,10)
polygon(c(0,x,qchisq(0.025,10)),c(0,y,0),col="lightblue")
points(x=qchisq(0.025,10),y=0,pch=19,cex=1)
x = seq(qchisq(0.975,10),30,0.01)
y = dchisq(x,10)
polygon(c(qchisq(0.975,10),x,30),c(0,y,0),col="lightblue")
points(x=qchisq(0.975,10),y=0,pch=19,cex=1)
text(9.8,0.02,expression(1-alpha),cex=1.25)
text(0.25,0.01,expression(alpha/2),cex=1.25)
text(25,0.0075,expression(alpha/2),cex=1.25)
text(15.5,0.08,expression(F(n[1]-1,n[2]-1)),cex=1.25)


# Example 7.3
# Page 138-141
rotary = c(127.75,127.87,127.86,127.92,128.03,127.94,127.91,128.10,128.01,128.11,127.79,127.93,
           127.89,127.96,127.80,127.94,128.02,127.82,128.11,127.92,127.74,127.78,127.85,127.96)
inline = c(127.90,127.90,127.74,127.93,127.62,127.76,127.63,127.93,127.86,127.73,127.82,127.84,
           128.06,127.88,127.85,127.60,128.02,128.05,127.95,127.89,127.82,127.92,127.71,127.78)
var(rotary) # sample variance
var(inline) # sample variance

# Figure 7.8
# Page 139
boxplot(rotary,inline,xlab="",names=c("Rotary","Inline"),ylab="Fill volume (in fluid ounces)",
	ylim=c(127.5,128.2),col="lightblue")

# F(23,23) pdf with 90% confidence
# Page 140
x = seq(0,4,0.001)
pdf = df(x,23,23)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",xaxt="n",yaxt="n",bty="n")
abline(h=0)
x = seq(0,qf(0.05,23,23),0.001)
y = df(x,23,23)
polygon(c(0,x,qf(0.05,23,23)),c(0,y,0),col="lightblue")
points(x=qf(0.05,23,23),y=0,pch=19,cex=1)
x = seq(qf(0.95,23,23),4,0.001)
y = df(x,23,23)
polygon(c(qf(0.95,23,23),x,4),c(0,y,0),col="lightblue")
points(x=qf(0.95,23,23),y=0,pch=19,cex=1)
text(1.7,0.75,expression(F(23,23)),cex=1.25)
text(1.15,0.2,0.95,cex=1.25)
text(0.05,0.1,0.05,cex=1.25)
text(2.7,0.08,0.05,cex=1.25)
text(qf(0.05,23,23),-0.0275,0.496,cex=1)
text(qf(0.95,23,23),-0.0275,2.014,cex=1)

# Confidence interval using var.test function in R
# Page 140
var.test(inline,rotary,conf.level=0.90)$conf.int

# Figure 7.9
# Page 141
# Left
library(car)
qqPlot(rotary,distribution="norm",mean=mean(rotary),sd=sd(rotary),
	xlab="Normal quantiles",ylab="Fill volume (in fluid ounces)",pch=16,
	envelope=list(border=TRUE,style="lines"),id=FALSE)
# Right
qqPlot(inline,distribution="norm",mean=mean(inline),sd=sd(inline),
	xlab="Normal quantiles",ylab="Fill volume (in fluid ounces)",pch=16,
	envelope=list(border=TRUE,style="lines"),id=FALSE)


# Figure 7.10
# Page 142
# Top--this is the same as Figure 7.7.
# Bottom left
x = seq(0,30,0.01)
pdf = dchisq(x,10)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",xaxt="n",yaxt="n",bty="n",ylim=c(0,0.1))
abline(h=0)
x = seq(0,qchisq(0.05,10),0.01)
y = dchisq(x,10)
polygon(c(0,x,qchisq(0.05,10)),c(0,y,0),col="lightblue")
points(x=qchisq(0.05,10),y=0,pch=19,cex=1)
text(9.8,0.02,expression(1-alpha),cex=1.25)
text(0.25,0.01,expression(alpha),cex=1.25)
text(15.5,0.08,expression(F(n[1]-1,n[2]-1)),cex=1.25)
# Bottom right
x = seq(0,30,0.01)
pdf = dchisq(x,10)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",xaxt="n",yaxt="n",bty="n",ylim=c(0,0.1))
abline(h=0)
x = seq(qchisq(0.95,10),30,0.01)
y = dchisq(x,10)
polygon(c(qchisq(0.95,10),x,30),c(0,y,0),col="lightblue")
points(x=qchisq(0.95,10),y=0,pch=19,cex=1)
text(9.8,0.02,expression(1-alpha),cex=1.25)
text(25,0.0075,expression(alpha),cex=1.25)
text(15.5,0.08,expression(F(n[1]-1,n[2]-1)),cex=1.25)


# Example 7.4
# Page 142-144
# Two-sided rejection region using alpha = 0.05
# Page 143
x = seq(0,4.5,0.001)
pdf = df(x,305,16)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",xaxt="n",yaxt="n",bty="n")
abline(h=0)
x = seq(0,qf(0.025,305,16),0.001)
y = df(x,305,16)
polygon(c(0,x,qf(0.025,305,16)),c(0,y,0),col="lightblue")
points(x=qf(0.025,305,16),y=0,pch=19,cex=1)
x = seq(qf(0.975,305,16),4,0.001)
y = df(x,305,16)
polygon(c(qf(0.975,305,16),x,4.5),c(0,y,0),col="lightblue")
points(x=qf(0.975,305,16),y=0,pch=19,cex=1)
text(1.8,0.8,expression(F(305,16)),cex=1.25)
text(1.2,0.2,0.95,cex=1.25)
text(0.06,0.1,0.025,cex=1.25)
text(2.9,0.08,0.025,cex=1.25)
text(qf(0.025,305,16),-0.0275,0.541,cex=1)
text(qf(0.975,305,16),-0.0275,2.342,cex=1)