# Tension Bond Strength data (Tab. 2-1, p. 26-28)
y1 <- c(16.85,16.40,17.21,16.35,16.52,17.04,16.96,17.15,16.59,16.57)
y2 <- c(16.62,16.75,17.37,17.12,16.98,16.87,17.34,17.02,17.08,17.27)
y <- data.frame(Modified=y1,Unmodified=y2)
y.means <- as.numeric(apply(y,2,mean))
opar <- par(mfrow=c(2,1),mar=c(5,7,4,2),las=1)
stripchart(y,xlab=expression("Strength (kgf/cm^2)"),pch=19)
arrows(y.means,rep(1.5,2),y.means,c(1.1,1.9),length=.1)
text(y.means,c(1.2,1.8),round(y.means,2),pos=4,cex=.8)

# Random deviates (instead of data from metal recovery used in the book)
rd <- rnorm(200,mean=70,sd=5)
hist(rd,xlab="quantile",ylab="Relative frequency",
main="Random Normal Deviates\n N(70,5)")
par(opar)

boxplot(y,ylab="Strength (kgf/cm^2)",las=1)

x <- seq(-3.5,3.5,by=0.01)
y <- dnorm(x)

plot(x,y,xlab="",ylab="",type="l",axes=F,lwd=2)
axis(side=1,cex.axis=.8); axis(2,pos=0,las=1,cex.axis=.8)
mtext(expression(mu),side=1,line=2,at=0)
mtext(expression(paste(frac(1, sigma*sqrt(2*pi)), " ",
plain(e)^{frac(-(x-mu)^2,
2*sigma^2)})),side=3,line=0)
# highlight a specific area (drawing is from left to right,
# then from right to left)
polygon(c(x[471:591],rev(x[471:591])),c(rep(0,121),rev(y[471:591])),
col="lightgray",border=NA)

# Table 2.2, p.41, two-sample t-test
# If we omit the var.equal=TRUE option, R computes the Welch modified t-test. 
t.test(y1,y2,var.equal=TRUE)
# c.f paired t test
t.test(y1,y2,paired=TRUE)

# one sided
t.test(y1,y2,var.equal=TRUE, alternative="less")
# c.f paired t test
t.test(y1,y2,paired=TRUE, alternative="less")
# one sided
t.test(y1,y2,var.equal=TRUE, alternative="greater")
# c.f paired t test
t.test(y1,y2,paired=TRUE, alternative="greater")

# calculate P value with given test statistics values e.g. problem 2.4, 6, 7
#for 2.4
z=c(2.35, 1.6, 2.15, -0.9, 1.645, -1.96)
twosidePvalue=2*(1-pnorm(abs(z)))
twosidePvalue
greaterAltPvalue=1-pnorm(z)
greaterAltPvalue
lessAltPvalue=pnorm(z)
lessAltPvalue

#for t test  use pt(t0, df) then.....