###############################################
## Author: Joshua M. Tebbs
## Date: 22 August 2025
## Update: 29 August 2025
## STAT 509 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
# See https://cran.r-project.org/web/packages/car/car.pdf for documentation

# Figure 7.1
# Page 107
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)


# Example 7.1
# Data analysis
# Page 108-109
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
t.test(bricks,conf.level=0.95)$conf.int


# Page 109
# t(35) pdf with quantiles
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,0.95,cex=1.25)
text(-3.5,0.04,0.025,cex=1.25)
text(3.5,0.04,0.025,cex=1.25)
text(1.5,0.3,"t(35)",cex=1.25)
text(2.21,-0.011,2.03,cex=1)


# Figure 7.2
# Page 112
# Left (Histogram)
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)
bins = seq(4,5,0.1)
hist(bricks,breaks=bins,xlab="Weight (in lbs)",ylab="Count",ylim=c(0,15),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 7.3
# Page 113
q = seq(0,45,0.01)
plot(q,dchisq(q,5),type="l",lty=1,xlab="q",ylab=expression(f[Q](q)),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.4
# Page 114
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)


# Example 7.2
# Data analysis
# Page 115-117
screws = c(1.194,1.177,1.204,1.195,1.209,1.208,1.210,1.206,1.187,1.187,1.219,1.198,1.196,1.194,1.194,
	1.208,1.207,1.201,1.214,1.203,1.198,1.195,1.207,1.185,1.189,1.191,1.204,1.199,1.196,1.212,
	1.198,1.188,1.203,1.199,1.211,1.215,1.202,1.206,1.212,1.189)
options(digits=5)
var(screws) # sample variance
# CI for population variance function
var.ci = function(data,conf.level=0.99){
    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(screws,conf.level=0.99) # CI for population variance
options(digits=1)
sqrt(var.ci(screws,conf.level=0.99)) # CI for population standard deviation


# Page 109
# chi^2(39) pdf with quantiles
x = seq(0,30,0.001)
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.005,10),0.001)
y = dchisq(x,10)
polygon(c(0,x,qchisq(0.005,10)),c(0,y,0),col="lightblue")
points(x=qchisq(0.005,10),y=0,pch=19,cex=1)
x = seq(qchisq(0.995,10),30,0.001)
y = dchisq(x,10)
polygon(c(qchisq(0.995,10),x,30),c(0,y,0),col="lightblue")
points(x=qchisq(0.995,10),y=0,pch=19,cex=1)
text(9.8,0.025,0.99,cex=1.25)
text(0.18,0.01,0.005,cex=1.25)
text(28,0.005,0.005,cex=1.25)
text(13.5,0.075,expression(paste(chi^2, "(39)")),cex=1.25)
text(2,-0.002,19.996,cex=1)
text(25,-0.002,65.476,cex=1)


# Figure 7.5
# Page 118
screws = c(1.194,1.177,1.204,1.195,1.209,1.208,1.210,1.206,1.187,1.187,1.219,1.198,1.196,1.194,1.194,
	1.208,1.207,1.201,1.214,1.203,1.198,1.195,1.207,1.185,1.189,1.191,1.204,1.199,1.196,1.212,
	1.198,1.188,1.203,1.199,1.211,1.215,1.202,1.206,1.212,1.189)
# Left (Histogram)
bins = seq(1.8,2.2,0.2)
hist(screws,breaks=bins,xlab="Diameter (in cm)",ylab="Count",ylim=c(0,15),main="",col="lightblue")
# Right (QQ plot)
qqPlot(screws,distribution="norm",mean=mean(screws),sd=sd(screws),
	xlab="Normal quantiles",ylab="Diameter (in cm)",pch=16,
	envelope=list(border=TRUE,style="lines"),id=FALSE)


# Figure 7.6
# Page 121
# Left
n = 100 # sample size 
B = 10000 # number of samples we will simulate
p = 0.30 # population proportion 
binomial.data = rbinom(B,n,p) 
sample.prop = binomial.data/n    
bins = seq(0.10,0.50,0.02)  
hist(sample.prop,prob=TRUE,breaks=bins,xlab="Sample proportion",ylab="Density",main = "",col="lightblue",xaxp=c(0.10,0.50,8),ylim=c(0,10))
# Add estimate of population density curve 
lines(sort(sample.prop),dnorm(sort(sample.prop),0.3,sqrt(0.3*0.7/100)),col="red",lwd=2)
# Right
n = 1000 # sample size 
B = 10000 
p = 0.30
binomial.data = rbinom(B,n,p) 
sample.prop = binomial.data/n    
bins = seq(0.24,0.36,0.005)  
hist(sample.prop,prob=TRUE,breaks=bins,xlab="Sample proportion",ylab="Density",main = "",col="lightblue",xaxp=c(0.24,0.36,12),ylim=c(0,30))
lines(sort(sample.prop),dnorm(sort(sample.prop),0.3,sqrt(0.3*0.7/1000)),col="red",lwd=2)


# Figure 7.7
# Page 122
# This figure is the same as Figure 7.1.


# Page 124
# N(0,1) pdf with quantiles
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,0.95,cex=1.25)
text(-3.5,0.04,0.025,cex=1.25)
text(3.5,0.04,0.025,cex=1.25)
text(1.5,0.3,"N(0,1)",cex=1.25)
text(2.21,-0.011,1.96,cex=1)


# Example 7.3
# Data analysis
# Page 124-125
# CI for population proportion function
p.ci = function(x,n,conf.level=0.95){
      est = x/n
      se = sqrt(est*(1-est)/n)
      z.upper = qnorm((1+conf.level)/2,0,1)
      c(est-z.upper*se,est+z.upper*se)
      }
options(digits=2)
p.ci(10,74,conf.level=0.95)