###############################################
## Author: Joshua M. Tebbs
## Date: 27 July 2016
## Update: 15 March 2024
## STAT 110 course notes: R Code Chapter 13
###############################################

# Example 13.1
# Page 119
# Turtle data
# Enter data
shell = c(26.96,33.42,26.57,23.50,36.46,40.54,36.16,33.32,35.99,
               43.55,44.33,42.73,42.15,47.43,49.96,46.04,48.76,47.78,
               45.81,49.05,17.65,44.71,54.29,52.01,51.15,54.42,52.62,
               53.27,54.07,50.40,59.69,56.30,54.29,54.58,55.11,57.65,
               56.35,55.68,58.40,58.06,57.79,56.54,57.03,57.64,59.27,
               64.79,61.96,60.08,9.34,73.84,60.61,64.91,60.35,62.63,
               63.33,63.00,64.55,66.03,64.75,60.24,69.01,65.07,65.77,
               65.30,68.54,65.28,67.54,68.49,66.98,65.67,70.26,70.94,
               73.52,71.01,77.34,74.63)
# Make histogram
# Figure 13.1 (Page 119)
bins = seq(5,85,5)
hist(shell,breaks=bins,xlab="Curved shell length (in cm)",ylab="Count",main="",col="lightblue",xlim=c(0,90))

# Make histogram with density curve superimposed
# Figure 13.2 (Page 121)
bins = seq(5,85,5)
hist(shell,breaks=bins,prob=TRUE,xlab="Curved shell length (in cm)",ylab="Density",main="",col="lightblue",xlim=c(0,90))
lines(density(shell.length),col="red",lwd=2)


# Figure 13.3
# Page 122
# BMI population density curve
y = seq(5,40,0.1)
pdf = dgamma(y,21.83,scale=0.92)
plot(y,pdf,type="l",xlab="BMI",ylab="Population density curve")
abline(h=0)


# Figure 13.4
# Page 123
# Each subfigure is created separately

# Upper left
x = seq(5,40,0.1)
pdf = dgamma(x,21.83,scale=0.92)
plot(x,pdf,type="l",xlab="BMI",main = "Underweight",ylab="Population density curve")
abline(h=0)
x = seq(5,14.5,0.001)
y = dgamma(x,21.83,scale=0.92)
polygon(c(0,x,14.5),c(0,y,0),col="lightblue")
points(x=14.5,y=0,pch=19,cex=1.5)
text(9.8,0.02,0.085,cex = 1.25)
# Calculate proportion (area)
pgamma(14.5,21.83,scale=0.92)

# Upper right
x = seq(5,40,0.1)
pdf = dgamma(x,21.83,scale=0.92)
plot(x,pdf,type="l",xlab="BMI",main = "Healthy",ylab="Population density curve")
abline(h=0)
x = seq(14.5,19.5,0.001)
y = dgamma(x,21.83,scale=0.92)
polygon(c(14.5,x,19.5),c(0,y,0),col="lightblue")
points(x=14.5,y=0,pch=19,cex=1.5)
points(x=19.5,y=0,pch=19,cex=1.5)
text(12.8,0.07,0.389,cex = 1.25)
# Calculate proportion (area)
pgamma(19.5,21.83,scale=0.92)-pgamma(14.5,21.83,scale=0.92)

# Lower left
x = seq(5,40,0.1)
pdf = dgamma(x,21.83,scale=0.92)
plot(x,pdf,type="l",xlab="BMI",main = "Overweight",ylab="Population density curve")
abline(h=0)
x = seq(19.5,21.5,0.001)
y = dgamma(x,21.83,scale=0.92)
polygon(c(19.5,x,21.5),c(0,y,0),col="lightblue")
points(x=19.5,y=0,pch=19,cex=1.5)
points(x=21.5,y=0,pch=19,cex=1.5)
text(24,0.09,0.179,cex = 1.25)
# Calculate proportion (area)
pgamma(21.5,21.83,scale=0.92)-pgamma(19.5,21.83,scale=0.92)

# Lower right
x = seq(5,40,0.1)
pdf = dgamma(x,21.83,scale=0.92)
plot(x,pdf,type="l",xlab="BMI",main = "Obese",ylab="Population density curve")
abline(h=0)
x = seq(21.5,40,0.001)
y = dgamma(x,21.83,scale=0.92)
polygon(c(21.5,x,40),c(0,y,0),col="lightblue")
points(x=21.5,y=0,pch=19,cex=1.5)
text(28.5,0.04,0.347,cex = 1.25)
# Calculate proportion (area)
1-pgamma(21.5,21.83,scale=0.92)


# Figure 13.5
# Page 124
# Each subfigure is created separately

# Upper left
y = seq(0,3,0.01)
pdf = dnorm(y,1.5,sqrt(0.16))
plot(y,pdf,type="l",xlab="",ylab="")
abline(h=0)

# Upper right
y = seq(0,80,0.1)
pdf = dgamma(y,4,1/6)
plot(y,pdf,type="l",xlab="",ylab="")
abline(h=0)

# Lower left
y = seq(0,1,0.01)
pdf = dbeta(y,5,2)
plot(y,pdf,type="l",xlab="",ylab="")
abline(h=0)

# Lower right
y = seq(0,110,0.1)
pdf = 0.25*dnorm(y,30,10) + 0.75*dnorm(y,70,12)
plot(y,pdf,type="l",xlab="",ylab="")
abline(h=0)


# Figure 13.6
# Page 128
# Left
x = seq(150,850,1)
pdf = dnorm(x,500,100)
plot(x,pdf,type="l",xlab="SAT math score",ylab="",xaxp=c(200,800,6))
points(x=500,y=0,pch=19,cex=1.5)
points(x=650,y=0,pch=4,cex=1.5)
abline(h=0)
# Right
x = seq(-2,38,0.11)
pdf = dnorm(x,18,6)
plot(x,pdf,type="l",xlab="ACT math score",ylab="",xaxp=c(0,36,6))
points(x=18,y=0,pch=19,cex=1.5)
points(x=21,y=0,pch=4,cex=1.5)
abline(h=0)


# Figure 13.8
# Page 130
x = seq(70,180,0.1)
pdf = dnorm(x,125,15)
plot(x,pdf,type="l",xlab="Systolic blood pressure (mm Hg)",ylab="",xaxp=c(80,170,6))
abline(h=0)


# Figure 13.9
# Page 133
x = seq(70,180,0.1)
pdf = dnorm(x,125,15)
plot(x,pdf,type="l",xlab="Systolic blood pressure (mm Hg)",ylab="",xaxp=c(80,170,6))
abline(h=0)
# Shading
x = seq(70,100,0.1)
y = dnorm(x,125,15)
polygon(c(70,x,100),c(0,y,0),col="lightblue")
points(x=100,y=0,pch=19,cex=1.5)
text(80,0.004,0.0446,cex=1.25)


# Figure 13.10
# Page 135
x = seq(70,180,0.1)
pdf = dnorm(x,125,15)
plot(x,pdf,type="l",xlab="Systolic blood pressure (mm Hg)",ylab="",xaxp=c(80,170,6))
abline(h=0)
# Shading
x = seq(130,180,0.1)
y = dnorm(x,125,15)
polygon(c(130,x,180),c(0,y,0),col="lightblue")
points(x=130,y=0,pch=19,cex=1.5)
text(160,0.01,0.3821,cex=1.25)


# Figure 13.11
# Page 136
x = seq(70,180,0.1)
pdf = dnorm(x,125,15)
plot(x,pdf,type="l",xlab="Systolic blood pressure (mm Hg)",ylab="",xaxp=c(80,170,6))
abline(h=0)
# Shading
x = seq(70,100,0.1)
y = dnorm(x,125,15)
polygon(c(70,x,100),c(0,y,0),col="lightblue")
points(x=100,y=0,pch=19,cex=1.5)
text(80,0.004,0.0446,cex=1.25)
x = seq(130,180,0.1)
y = dnorm(x,125,15)
polygon(c(130,x,180),c(0,y,0),col="lightblue")
points(x=130,y=0,pch=19,cex=1.5)
text(160,0.01,0.3821,cex=1.25)


# Figure 13.12
# Page 137
x = seq(115,385,0.1)
pdf = dnorm(x,250,40)
plot(x,pdf,type="l",xlab="Marathon completion time (in minutes)",ylab="",xaxp=c(130,370,6))
abline(h=0)