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

# Example 6.1
# Battery lifetime data
# Page 69-70
battery = c(4285,2066,2584,1009,318,1429,981,1402,1137,414, 
	564,604,14,4152,737,852,1560,1786,520,396,
	1278,209,349,478,3032,1461,701,1406,261,83,
	205,602,3770,726,3894,2662,497,35,2778,1379, 
	3920,1379,99,510,582,308,3367,99,373,454)
# Page 70
# Create histogram
hist(battery,xlab="Lifetime (in hours)",xlim=c(0,5000),ylab="Count",main="",col="lightblue",cex.lab=1.25)
# Create boxplot
# Use range=1.5 to detect outliers
boxplot(battery,range=0,xlab="",ylab="Lifetime (in hours)",col="lightblue",cex.lab=1.25)
# Calculate summary statistics
# Page 72
mean(battery) ## sample mean
var(battery) ## sample variance
sd(battery) ## sample standard deviation


# Example 6.2
# Brake time example
# Page 75-76
y = seq(0,3,0.01)
# Plot population distribution
plot(y,dnorm(y,1.5,sqrt(0.16)),type="l",lty=1,xlab="Brake times (in seconds)",ylab="PDF",ylim=c(0,5))
# Add sampling distribution when n=5
lines(y,dnorm(y,1.5,sqrt(0.16/5)),lty=4)
# Add sampling distribution when n=25
lines(y,dnorm(y,1.5,sqrt(0.16/25)),lty=8)
abline(h=0)
# Add legend
legend(0,5,lty = c(1,4,8), 
	c(
	expression(paste("Population distribution")),
	expression(paste("Sample mean, n=5")),
	expression(paste("Sample mean, n=25"))
	))


# Example 6.3
# Rat example
# Page 77-78
y = seq(0,20,0.01)
# Plot population distribution
plot(y,dexp(y,1/5),type="l",lty=1,xlab="Time to death (in days)",ylab="PDF",ylim=c(0,0.4))
# Add sampling distribution when n=5
lines(y,dgamma(y,5,1),lty=4)
# Add sampling distribution when n=5
lines(y,dgamma(y,25,5),lty=8)
abline(h=0)
abline(v=0,lty=2)
# Add legend
legend(10,0.3,lty = c(1,4,8), 
	c(
	expression(paste("Population distribution")),
	expression(paste("Sample mean, n=5")),
	expression(paste("Sample mean, n=25"))
	))


# Figure 6.4
# Page 80
# N(0,1) and t PDFs
y = seq(-5,5,0.01)
# Plot N(0,1) pdf
plot(y,dnorm(y,0,1),type="l",lty=1,xlab="",ylab="PDF",ylim=c(0,0.4))
# Add other pdfs
lines(y,dt(y,2),lty=4)
lines(y,dt(y,10),lty=8)
abline(h=0)
# Add legend
legend(2,0.35,lty = c(1,4,8), 
	c(
	expression(paste("N(0,1)")),
	expression(paste("t(2)")),
	expression(paste("t(10)"))
	))


# Example 6.5
# Page 81-84
# Pipe diameter data
pipes = c(1.296,1.320,1.311,1.298,1.315,1.305,1.278,1.294,1.311,1.290,1.284,1.287,1.289,1.292,1.301,
	1.298,1.287,1.302,1.304,1.301,1.313,1.315,1.306,1.289,1.291)
# Calculate summary statistics
mean(pipes) ## sample mean
sd(pipes) ## sample standard deviation

# Page 82
# Plot t(24) pdf
t = seq(-5,5,0.001)
pdf = dt(t,24)
plot(t,pdf,type="l",lty=1,xlab="t",ylab="PDF")
abline(h=0)
# Add points
points(x=4.096,y=0,pch=4,cex=1.5)

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