###############################################
## Author: Joshua M. Tebbs
## Date: 5 July 2026
## Update: 16 July 2026
## STAT 515 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 

# Install and load "car" R package
# For normal qq plots
# See https://cran.r-project.org/web/packages/car/car.pdf for documentation

# Figure 4.1
# Page 3
bottles = rbinom(5000,1,0.10)
# Create plot of the proportion of defective bottles
# Initialize
prop = bottles*0
for (i in 1:5000){
               prop[i] = sum(bottles[1:i])/i
               }
plot(prop,ylab="Proportion of defective bottles",xlab="Number of bottles tested",type="o",pch=16)
abline(h=0.10,lty=2)
prop[5000] # estimate of P(E) after 5000 simulated bottles (this will change from simulation to simulation)


# Figure 4.2
# Page 15
y = c(1,2,3,4,5)
prob = c(0.38,0.27,0.18,0.11,0.06)
plot(y,prob,type="h",xlab="y",ylab=expression(p[Y](y)),ylim=c(0,0.4),cex.lab=1.25,lwd=2)
abline(h=0)
# Simulation
options(digits=3) # control number of significant digits presented
n = 500
y = c(1,2,3,4,5)
prob = c(0.38,0.27,0.18,0.11,0.06)
rating = sample(y,n,replace=TRUE,prob=prob)
mean(rating)


# Figure 4.3
# Page 19
y = c(0,1,2,3,4)
prob = dbinom(y,4,0.4)
plot(y,prob,type="h",xlab="y",ylab=expression(p[Y](y)),xaxp=c(0,4,4),ylim=c(0,0.4),cex.lab=1.25,lwd=2)
abline(h=0)


# Figure 4.4
# Page 21
y = seq(0,30,1)
prob = dbinom(y,30,0.05)
plot(y,prob,type="h",xlab="y",ylab=expression(p[Y](y)),ylim=c(0,max(prob)),cex.lab=1.25,lwd=2)
abline(h=0)


# Figure 4.5
# Page 23
y = seq(0,10,1)
prob = dpois(y,1.5)
plot(y,prob,type="h",xlab="y",ylab=expression(p[Y](y)),ylim=c(0,max(prob)),cex.lab=1.25,lwd=2)
abline(h=0)


# Figure 4.6
# Page 26
y = seq(0,30,0.01)
pdf = 3000/(10+y)^4
plot(y,pdf,type="l",xlab="y",ylab=expression(f[Y](y)),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.7
# Page 27
y = seq(4.9,5.1,0.01)
pdf = 5+y-y
plot(y,pdf,type="l",xlab="y",ylab=expression(f[Y](y)),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.8
# Page 28
y = seq(-10,10,0.01)
plot(y,dnorm(y,0,1),type="l",lty=1,xlab="y",ylab=expression(f[Y](y)),cex.lab=1.25)
lines(y,dnorm(y,-2,2),lty=4)
lines(y,dnorm(y,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.9
# Page 29
y = seq(350,550,0.05)
pdf = dnorm(y,450,30)
plot(y,pdf,type="l",xlab="y",ylab=expression(f[Y](y)),xaxp=c(360,540,6),cex.lab=1.25)
abline(h=0)


# Figure 4.10
# Page 31
y = seq(75,175,0.05)
pdf = dnorm(y,125,15)
plot(y,pdf,type="l",xlab="y",ylab=expression(f[Y](y)),xaxp=c(80,170,6),cex.lab=1.25)
abline(h=0)


# Figure 4.11
# Page 33
y = seq(-3.5,3.5,0.001)
pdf = dnorm(y)
plot(y,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.2),0.001)
y = dnorm(x)
polygon(c(-3.5,x,qnorm(0.2)),c(0,y,0),col="lightblue")
points(x=qnorm(0.2),y=0,pch=19,cex=1)
text(-1.5,0.05,expression(p),cex=1.25)
text(0.5,0.1,expression(1-p),cex=1.25)


# Figure 4.12
# Page 34
y = seq(115,385,0.1)
pdf = dnorm(y,250,40)
plot(y,pdf,type="l",xlab="y",ylab=expression(f[Y](y)),xaxp=c(130,370,6),cex.lab=1.25)
abline(h=0)


# Figure 4.13
# Page 36
# Upper left
pennies = rdunif(1000,0,40) # generate population of penny ages
bins = seq(0,40,1) 
hist(pennies,breaks=bins,xlab="Age (in years)",ylab="Count",main="",col="lightblue") # population distribution

# Lower right
# Sampling distribution of sample mean with n=20
n = 20 # sample size
B = 100000 # number of samples
storage = matrix(0,nrow=B,ncol=n) # initialize
# This for loop creates B samples of size n
for (i in 1:B){
	storage[i,] = sample(pennies,n,replace=F)
	}
sample.mean = rep(0,B) # initialize
for (j in 1:B){
	sample.mean[j] = mean(storage[j,])
	}
sample.mean 
hist(sample.mean,breaks=bins,xlab="Sample mean (n=20)",ylab="Count",main="",col="lightblue")


# Figure 4.14
# Page 37
# Lower right
# Sampling distribution of sample variance with n=20
n = 20 # sample size
B = 100000 # number of samples
storage = matrix(0,nrow=B,ncol=n) # initialize
# This for loop creates B samples of size n
for (i in 1:B){
	storage[i,] = sample(pennies,n,replace=F)
	}
sample.var = rep(0,B) # initialize
for (j in 1:B){
	sample.var[j] = var(storage[j,])
	} 
hist(sample.var,xlab="Sample variance (n=20)",ylab="Count",main="",col="lightblue")


# Figure 4.15
# Page 40
y = seq(28,84,0.1)
# Plot population distribution
plot(y,dnorm(y,56,8),type="l",lty=1,xlab="Ejection fraction",ylab="",xaxp=c(32,80,6),ylim=c(0,0.2))
# Add sampling distribution
lines(y,dnorm(y,56,2),lty=4)
abline(h=0)
points(x=56,y=0,pch=19,cex=1.25)
# Add legend
legend(27.5,0.2,lty = c(1,4), 
	c(
	expression(paste("Population distribution")),
	expression(paste("Sampling distribution"))
	))


# Figure 4.16
# Page 42
# Left (Poisson PMF)
y = seq(0,10,1)
prob = dpois(y,1.5)
plot(y,prob,type="h",xlab="y",ylab=expression(p[Y](y)),ylim=c(0,max(prob)),cex.lab=1.25,lwd=2)
abline(h=0)
points(x=1.5,y=0,pch=19,cex=1.25)
# Right (Sampling distribution of sample mean)
x = seq(0.8,2.2,0.0001)
pdf = dnorm(x,1.5,sqrt(0.03))
plot(x,pdf,type="l",xlab="Sample mean number of defects (n = 50)",xaxp=c(0.75,2.25,6),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.0019,cex=1)


# Figure 4.17
# Page 46
y = seq(420,580,1)
pdf = dbinom(y,1000,0.5)
plot(y,pdf,type="h",xlab="y",ylab="Binomial pmf",xaxp=c(420,580,4),ylim=c(0,max(pdf)),cex.lab=1.25)
abline(h=0)
x = seq(400,600,0.01)
lines(x,dnorm(x,500,sqrt(250)),lty=1,col="red")


# Figure 4.18
# Page 48
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 4.19
# Page 49
library(car)
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 4.20
# Page 50
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)
library(car)
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)