###############################################
## Author: Joshua M. Tebbs
## Date: 13 January 2019
## Update: 18 February 2019
## STAT 512 course notes: R Code
## Chapter 8
###############################################

# Figure 8.1
# Page 75
# Left
y = seq(10,30,0.01)
plot(y,dnorm(y,20,2.5),type="l",lty=1,bty="n",xaxt="n",yaxt="n",xlab="",ylab="",ylim=c(-0.01,0.16),cex.lab=1.25)
abline(h=0)
points(x=20,y=0,pch=19,cex=1.25)
text(20,-0.01,expression(theta),cex=1.5)
text(30,-0.0075,expression(hat(theta)),cex=1.5)
# Right
y = seq(10,30,0.01)
plot(y,dnorm(y,20,2.5),type="l",lty=1,bty="n",xaxt="n",yaxt="n",xlab="",ylab="",ylim=c(-0.01,0.16),cex.lab=1.25)
abline(h=0)
points(x=17,y=0,pch=19,cex=1.25)
text(17,-0.01,expression(theta),cex=1.5)
text(30,-0.0075,expression(hat(theta)),cex=1.5)


# Figure 8.2
# Page 76
y = seq(0,1,0.01)
pdf = rep(1,length(y))
plot(y,pdf,type="l",xlab="y",ylab="PDF",xaxt='n',yaxt='n',xlim=c(0,1.25),ylim=c(0,1.25),cex.lab=1.25)
axis(side=1,at=c(0,1),labels=expression(0,theta))
axis(side=2,at=c(0,1),labels=expression(0,1/theta))
lines(c(1,1.5),c(0,0),lty=1)


# Figure 8.3
# Page 78
# Left
y = seq(10,30,0.01)
plot(y,dnorm(y,20,2.5),type="l",lty=1,bty="n",xaxt="n",yaxt="n",xlab="",ylab="",ylim=c(-0.01,0.16),cex.lab=1.25)
abline(h=0)
points(x=20,y=0,pch=19,cex=1.25)
text(20,-0.01,expression(theta),cex=1.5)
text(30,-0.0075,expression(hat(theta)[1]),cex=1.5)
# Right
y = seq(10,30,0.01)
plot(y,dnorm(y,20,0.5),type="l",lty=1,bty="n",xaxt="n",yaxt="n",xlab="",ylab="",ylim=c(-0.04,0.8),cex.lab=1.25)
abline(h=0)
points(x=20,y=0,pch=19,cex=1.25)
text(20,-0.04,expression(theta),cex=1.5)
text(30,-0.035,expression(hat(theta)[2]),cex=1.5)


# Figure 8.4
# Page 80
n = 25
p = seq(0,1,0.001)
mse.1 = p*(1-p)/n
mse.2 = (n*p*(1-p))/((n+4)^2) + ((n*p+2)/(n+4)-p)^2
plot(p,mse.1,type="l",lty=1,xlab="p",ylab="MSE",xlim=c(0,1),cex.lab=1.25)
abline(h=0)
lines(p,mse.2,lty=4)
legend(0.37,0.004,lty = c(1,4),c(expression(MSE(hat(p)[1])),expression(MSE(hat(p)[2]))))


# Figure 8.5
# Page 85
y = seq(10,30,0.01)
plot(y,dnorm(y,20,2.5),type="l",bty="n",lty=1,xaxt="n",yaxt="n",xlab="",ylab="",ylim=c(-0.01,0.16),cex.lab=1.25)
abline(h=0)
points(x=20,y=0,pch=19,cex=1.25)
text(20,-0.01,expression(mu),cex=1.5)
text(30,-0.0075,expression(bar(y)),cex=1.5)
x = seq(25,30,0.01)
y = dnorm(x,20,2.5)
polygon(c(25,x,30),c(0,y,0),col="lightblue")
x = seq(10,15,0.01)
y = dnorm(x,20,2.5)
polygon(c(10,x,15),c(0,y,0),col="lightblue")
text(20,0.07,0.95,cex=1.25)


# Figure 8.6
# Page 87
y = seq(10,30,0.01)
plot(y,dnorm(y,20,2.5),type="l",bty="n",lty=1,xaxt="n",yaxt="n",xlab="",ylab="",ylim=c(-0.01,0.16),cex.lab=1.25)
abline(h=0)
points(x=20,y=0,pch=19,cex=1.25)
text(20,-0.01,expression(p),cex=1.5)
text(30,-0.0075,expression(hat(p)),cex=1.5)
x = seq(25,30,0.01)
y = dnorm(x,20,2.5)
polygon(c(25,x,30),c(0,y,0),col="lightblue")
x = seq(10,15,0.01)
y = dnorm(x,20,2.5)
polygon(c(10,x,15),c(0,y,0),col="lightblue")
text(20,0.07,0.95,cex=1.25)


# Figure 8.7
# Page 91
birth.weights = read.table("C:\\Users\\tebbs.MATHSTAT\\Documents\\texfiles\\Classes\\USC\\stat512\\s19\\data\\birthweights.txt",header=TRUE)
birth.weights = birth.weights[,1]
hist(birth.weights,xlab="Birth weight (in grams)",ylab="Count",main="",col="lightblue",ylim=c(0,120))


# Figure 8.8
# Page 93
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.5)
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.5)
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)


# Figure 8.9
# Page 94
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.025,10),0.001)
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.5)
x = seq(qchisq(0.975,10),30,0.001)
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.5)
text(9.8,0.02,expression(1-alpha),cex=1.25)
text(0.5,0.01,expression(alpha/2),cex=1.25)
text(25,0.0075,expression(alpha/2),cex=1.25)


# Figure 8.10
# Page 97
x = seq(0,1,0.001)
pdf = dbeta(x,2.5,1)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",xaxt="n",yaxt="n",bty="n",xlim=c(0,1.05))
abline(h=0)
x = seq(0,qbeta(0.025,2.5,1),0.001)
y = dbeta(x,2.5,1)
polygon(c(0,x,qbeta(0.025,2.5,1)),c(0,y,0),col="lightblue")
points(x=qbeta(0.025,2.5,1),y=0,pch=19,cex=1.5)
x = seq(qbeta(0.975,2.5,1),1,0.001)
y = dbeta(x,2.5,1)
polygon(c(qbeta(0.975,2.5,1),x,1),c(0,y,0),col="lightblue")
points(x=qbeta(0.975,2.5,1),y=0,pch=19,cex=1.5)
text(0.7,0.6,expression(1-alpha),cex=1.25)
text(0.1,0.25,expression(alpha/2),cex=1.25)
text(1.05,0.25,expression(alpha/2),cex=1.25)


# Figure 8.11
# Page 103
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.5)
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.5)
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)


# Figure 8.12
# Page 105
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.025,10),0.001)
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.5)
x = seq(qchisq(0.975,10),30,0.001)
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.5)
text(9.8,0.02,expression(1-alpha),cex=1.25)
text(0.5,0.01,expression(alpha/2),cex=1.25)
text(25,0.0075,expression(alpha/2),cex=1.25)


# Figure 8.13
# 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.5)
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.5)
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 8.15
# Pages 108-109
loc.1 = c(5030,4980,13700,11910,10730,8130,11400,26850,860,17660,2200,22800,4250,1130,15040,1690)
loc.2 = c(2800,2810,4670,1330,6890,3320,7720,1230,7030,2130,7330,2190)
t.test(loc.1,loc.2,conf.level=0.95,var.equal=TRUE)$conf.int
t.test(loc.1,loc.2,conf.level=0.95,var.equal=FALSE)$conf.int


# Figure 8.14
# Page 110
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.025,10),0.001)
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.5)
x = seq(qchisq(0.975,10),30,0.001)
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.5)
text(9.8,0.02,expression(1-alpha),cex=1.25)
text(0.5,0.01,expression(alpha/2),cex=1.25)
text(25,0.0075,expression(alpha/2),cex=1.25)