###############################################
## Author: Joshua M. Tebbs
## Date: 11 December 2011
## Update: 30 December 2017
## STAT 509 course notes: R Code
## Chapter 7
###############################################

# Figure 7.1
# Page 87
# t pdf with quantiles and shading
x = seq(-4,4,0.001)
pdf = dt(x,10)
plot(x,pdf,type="l",lty=1,xlab="t",ylab="PDF",ylim=c(0,0.4))
abline(h=0)
x = seq(-4,qt(0.025,10),0.001)
y = dt(x,10)
polygon(c(-4,x,qt(0.025,10)),c(0,y,0),col="grey")
# Add points
points(x=qt(0.025,10),y=0,pch=19,cex=1.5)
x = seq(qt(0.975,10),4,0.001)
y = dt(x,10)
polygon(c(qt(0.975,10),x,4),c(0,y,0),col="grey")
# Add points
points(x=qt(0.975,10),y=0,pch=19,cex=1.5)
# Add text
text(-0.025,0.1,expression(1-alpha))
text(-3.5,0.05,expression(alpha/2))
text(3.5,0.05,expression(alpha/2))


# Example 7.1
# Cadmium data
# Page 89-92
cadmium = c(0.044,0.030,0.052,0.044,0.046,0.020,0.066,
	0.052,0.049,0.030,0.040,0.045,0.039,0.039,	
	0.039,0.057,0.050,0.056,0.061,0.042,0.055,
	0.037,0.062,0.062,0.070,0.061,0.061,0.058,
	0.053,0.060,0.047,0.051,0.054,0.042,0.051)
# Calculate summary statistics
mean(cadmium) ## sample mean
sd(cadmium) ## sample standard deviation

# Calculate confidence interval directly
# Page 91
t.test(cadmium,conf.level=0.99)$conf.int

# Create normal qqplot
# Page 92
# Need to install (then load) car package in R 
qqPlot(cadmium,distribution="norm",xlab="N(0,1) quantiles",ylab="Cadmium data",pch=16)


# Figure 7.3
# Page 94
# Chi-square PDFs
y = seq(0,45,0.01)
# Plot chi-square pdf with 5 df
plot(y,dchisq(y,5),type="l",lty=1,xlab=expression(chi^2),ylab="PDF",ylim=c(0,0.16))
# Add other pdfs
lines(y,dchisq(y,10),lty=4)
lines(y,dchisq(y,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 95
# chi-square pdf with quantiles and shading
x = seq(0,30,0.001)
pdf = dchisq(x,10)
plot(x,pdf,type="l",lty=1,xlab=expression(chi^2),ylab="PDF",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="grey")
# Add points
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="grey")
# Add points
points(x=qchisq(0.975,10),y=0,pch=19,cex=1.5)
# Add text
text(9.8,0.02,expression(1-alpha))
text(0.5,0.01,expression(alpha/2))
text(25,0.0075,expression(alpha/2))


# Example 7.2
# IKEA diameter data
# Page 97-99
diameters = c(1.206,1.190,1.200,1.195,1.201,1.200,1.198,1.196,1.195,
	1.202,1.203,1.210,1.206,1.193,1.207,1.201,1.199,1.200,1.199,
	1.204,1.194,1.203,1.194,1.199,1.203,1.200,1.197,1.208,1.199,
	1.205,1.199,1.204,1.202,1.196,1.211,1.204)

# R does not have an internal function to calculate CI for population variance (not that I could find quickly)
# I wrote this function to do it
var.interval = function(data,conf.level=0.95){
	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)
	}

# CI for population variance
var.interval(diameters)
# CI for population standard deviation
sd.interval = sqrt(var.interval(diameters))
sd.interval

# Create normal qqplot
# Page 99
# Need to install (then load) car package in R 
qqPlot(diameters,distribution="norm",xlab="N(0,1) quantiles",ylab="Diameter data",pch=16)


# Figure 7.6
# Page 102
# N(0,1) pdf with quantiles and shading
x = seq(-4,4,0.001)
pdf = dnorm(x,0,1)
plot(x,pdf,type="l",lty=1,xlab="z",ylab="PDF",ylim=c(0,0.4))
abline(h=0)
x = seq(-4,qnorm(0.025,0,1),0.001)
y = dnorm(x,0,1)
polygon(c(-4,x,qnorm(0.025,0,1)),c(0,y,0),col="grey")
# Add points
points(x=qnorm(0.025,0,1),y=0,pch=19,cex=1.5)
x = seq(qnorm(0.975,0,1),4,0.001)
y = dnorm(x,0,1)
polygon(c(qnorm(0.975,0,1),x,4),c(0,y,0),col="grey")
# Add points
points(x=qnorm(0.975,0,1),y=0,pch=19,cex=1.5)
# Add text
text(-0.025,0.1,expression(1-alpha))
text(-2.8,0.04,expression(alpha/2))
text(2.8,0.04,expression(alpha/2))