###############################################
## Author: Joshua M. Tebbs
## Date: 21 August 2018
## Update: 10 November 2018
## STAT 511 course notes: R Code
## Chapter 4
###############################################

# Figure 4.1
# Page 74
# Left
y = c(0,1,2,3)
prob = dbinom(y,3,0.4)
plot(y,prob,type="h",xlab="y",ylab="PMF",ylim=c(0,max(prob)),xaxt='n',cex.lab=1.25)
abline(h=0)
axis(side=1,at=y,labels=c("0","1","2","3"))
# Right
y = c(0,1,2,3)
cdf = c(0,cumsum(prob))
cdf.plot = stepfun(y,cdf,f=0)
plot.stepfun(cdf.plot,xlab="y",ylab="CDF",verticals=FALSE,do.points=TRUE,main="",pch=16,cex.lab=1.25)


# Figure 4.2
# Page 77
y = seq(0,4,0.01)
cdf = 1-exp(-y^2)
plot(y,cdf,type="l",xlab="y",ylab="CDF",xlim=c(-1,4),ylim=c(0,1),cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)


# Figure 4.3
# Page 78
bmi = read.table("C:\\Users\\Tebbs\\Documents\\texfiles\\Classes\\USC\\stat511\\f18\\data\\bmi.txt",header=TRUE)
bmi = bmi[,1]
bins = seq(10,50,2)
hist(bmi,prob=TRUE,breaks=bins,xlab="Body mass index (BMI)",ylab="Proportion",main="",col="lightblue",cex.lab=1.25)
lines(density(bmi),col="red",lwd=2)


# Figure 4.4
# Page 80
# Left
y = seq(0,2,0.01)
pdf = (3/8)*(4*y-2*y^2)
plot(y,pdf,type="l",xlab="y",ylab="PDF",xlim=c(-1,3),cex.lab=1.25)
abline(h=0)
# Right
y = seq(0,2,0.01)
cdf = c((3/8)*(2*y^2 - (2/3)*y^3))
plot(y,cdf,type="l",xlab="y",ylab="CDF",xlim=c(-1,3),ylim=c(0,1),cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)
lines(c(2,3),c(1,1),lty=1)


# Figure 4.5
# Page 81
# Left
y = seq(0,2,0.01)
pdf = (3/8)*(4*y-2*y^2)
plot(y,pdf,type="l",xlab="y",ylab="PDF",xlim=c(-1,3),cex.lab=1.25)
abline(h=0)
x = seq(0,0.5,0.001)
y = (3/8)*(4*x-2*x^2)
polygon(c(0,x,0.5),c(0,y,0),col="lightblue")
points(x=0.5,y=0,pch=19,cex=1.5)
# Right
y = seq(0,2,0.01)
cdf = c((3/8)*(2*y^2 - (2/3)*y^3))
plot(y,cdf,type="l",xlab="y",ylab="CDF",xlim=c(-1,3),ylim=c(0,1),cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)
lines(c(2,3),c(1,1),lty=1)
points(x=0.5,y=(3/8)*(2*(0.5)^2-(2/3)*(0.5)^3),pch=19,cex=1.5)


# Figure 4.6
# Page 82
# Left
y = seq(-7,7,0.01)
pdf = exp(-y)/(1+exp(-y))^2
plot(y,pdf,type="l",xlab="y",ylab="PDF",xlim=c(-7,7),cex.lab=1.25)
abline(h=0)
x = seq(-2,2,0.001)
y = exp(-x)/(1+exp(-x))^2
polygon(c(-2,x,2),c(0,y,0),col="lightblue")
points(x=-2,y=0,pch=19,cex=1.5)
points(x=2,y=0,pch=19,cex=1.5)
# Right
y = seq(-7,7,0.01)
cdf = 1/(1+exp(-y))
plot(y,cdf,type="l",xlab="y",ylab="CDF",ylim=c(0,1),cex.lab=1.25)
abline(h=0)
points(x=-2,y=1/(1+exp(2)),pch=19,cex=1.5)
points(x=2,y=1/(1+exp(-2)),pch=19,cex=1.5)


# Figure 4.7
# Page 84
# Left
y = seq(0,35,0.01)
pdf = 3000/(10+y)^4
plot(y,pdf,type="l",xlab="y",ylab="PDF",cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)
x = seq(0,10/((0.5)^(1/3))-10,0.001)
y = 3000/(10+x)^4
polygon(c(0,x,10/((0.5)^(1/3))-10),c(0,y,0),col="lightblue")
points(x=10/((0.5)^(1/3))-10,y=0,pch=19,cex=1.5)
# Right
y = seq(0,35,0.01)
cdf = 1-(10/(10+y))^3
plot(y,cdf,type="l",xlab="y",ylab="CDF",ylim=c(0,1),cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)
abline(h=0.5,lty=2)
points(x=10/((0.5)^(1/3))-10,y=0,pch=19,cex=1.5)


# Figure 4.8
# Page 86
y = seq(0,10,0.01)
pdf = y*exp(-y)
plot(y,pdf,type="l",xlab="y",ylab="PDF",cex.lab=1.25)
abline(h=0)
points(x=2,y=0,pch=19,cex=1.5)


# Figure 4.9
# Page 89
y = seq(0,6,0.01)
plot(y,dexp(y,1/0.5),type="l",lty=4,xlab="y",ylab="PDF",cex.lab=1.25)
lines(y,dexp(y,1),lty=1)
lines(y,dexp(y,1/2),lty=8)
abline(h=0)
abline(v=0,lty=2)
legend(3,1,lty = c(4,1,8),c(expression(paste(beta, "=0.5")),expression(paste(beta, "=1")),expression(paste(beta, "=2"))))


# Figure 4.10
# Page 91
# Left
y = seq(0,2*pi,0.01)
pdf = (1+y-y)/(2*pi)
matplot(y,pdf,xlab="y",ylab="PDF",type="l",xaxt="n",yaxt="n",xlim=c(0,7.5),ylim=c(0,0.5),cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)
axis(1, at = c(0, pi/2, pi, 3*pi/2, 2*pi),labels = expression(0, pi/2, pi, 3*pi/2, 2*pi))
axis(2, at = c(0, 1/(2*pi)),labels = expression(0, 1/2*pi))
# Right
cdf = y/(2*pi)
matplot(y,cdf,xlab="y",ylab="CDF",type="l",xaxt="n",xlim=c(0,7.5),ylim=c(0,1),cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)
lines(c(2*pi,7.5),c(1,1),lty=1)
axis(1, at = c(0, pi/2, pi, 3*pi/2, 2*pi),labels = expression(0, pi/2, pi, 3*pi/2, 2*pi))


# Figure 4.11
# Page 92
y = seq(-10,10,0.01)
plot(y,dnorm(y,0,1),type="l",lty=1,xlab="y",ylab="PDF",cex.lab=1.25)
lines(y,dnorm(y,-2,2),lty=4)
lines(y,dnorm(y,1,3),lty=8)
abline(h=0)
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.12
# Page 96
# Left
y = seq(70,180,0.1)
pdf = dnorm(y,125,15)
plot(y,pdf,type="l",xlab="y",ylab="PDF",xaxp=c(80,170,6),cex.lab=1.25)
abline(h=0)
x = seq(70,90,0.1)
y = dnorm(x,125,15)
polygon(c(70,x,90),c(0,y,0),col="lightblue")
points(x=90,y=0,pch=19,cex=1.5)
# Right
y = seq(70,180,0.1)
pdf = dnorm(y,125,15)
plot(y,pdf,type="l",xlab="y",ylab="PDF",xaxp=c(80,170,6),cex.lab=1.25)
abline(h=0)
points(x=137.6243,y=0,pch=19,cex=1.5)


# Figure 4.13
# Page 98
# Left
y = seq(0,1200,0.1)
pdf = dexp(y,1/190)
plot(y,pdf,type="l",xlab="y",ylab="PDF",cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)
x = seq(100,1200,0.1)
y = dexp(x,1/190)
polygon(c(100,x,1200),c(0,y,0),col="lightblue")
points(x=100,y=0,pch=19,cex=1.5)
# Right
y = seq(0,1200,0.1)
cdf = pexp(y,1/190)
plot(y,cdf,type="l",xlab="y",ylab="CDF",ylim=c(0,1),cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)
abline(h=0.9,lty=2)
points(x=437.5,y=0,pch=19,cex=1.5)


# Figure 4.14
# Page 101
y = seq(0,25,0.01)
plot(y,dgamma(y,1.5,1/1.6),type="l",lty=1,xlab="y",ylab="PDF",cex.lab=1.25)
lines(y,dgamma(y,2,1/2),lty=4)
lines(y,dgamma(y,2.5,1/3),lty=8)
abline(h=0)
legend(10,0.20,lty=c(1,4,8),c(expression(paste(alpha,"=1.5, ",beta,"=1.6")),expression(paste(alpha,"=2.0, ",beta,"=2.0")),expression(paste(alpha,"=2.5, ",beta,"=3.0"))))


# Figure 4.15
# Page 103
y = seq(0,1200,0.1)
pdf = dgamma(y,2.7,1/100)
plot(y,pdf,type="l",xlab="y",ylab="PDF",cex.lab=1.25)
abline(h=0)
x = seq(365,1200,0.1)
y = dgamma(x,2.7,1/100)
polygon(c(365,x,1200),c(0,y,0),col="lightblue")
points(x=365,y=0,pch=19,cex=1.5)


# Figure 4.16
# Page 106
y = seq(0,1,0.01)
plot(y,dbeta(y,3,3),type="l",lty=1,xlab="y",ylab="PDF",ylim=c(0,4),cex.lab=1.25)
lines(y,dbeta(y,2,6.5),lty=4)
lines(y,dbeta(y,4,1),lty=8)
abline(h=0)
abline(v=1,lty=2)
legend(0.4,3.5,lty=c(1,4,8),c(expression(paste(alpha,"=3, ",beta,"=3")),expression(paste(alpha,"=2, ",beta,"=6.5")),expression(paste(alpha,"=4, ",beta,"=1"))))


# Figure 4.17
# Page 108
y = seq(0,1,0.001)
pdf = dbeta(y,1,5)
plot(y,pdf,type="l",xlab="y",ylab="PDF",cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)
x = seq(qbeta(0.99,1,5),1,0.001)
y = dbeta(x,1,5)
polygon(c(qbeta(0.99,1,5),x,1),c(0,y,0),col="lightblue")
text(0.7,0.3,0.01,cex = 1)