# Chapter 8 Script for PASWR
# 11/09/07 ATA
#
library(PASWR)
# Example 8.1
# function interval.plot() is part of PASWR package
#
set.seed(402)
m<-100 # Number of samples
n<-500 # Sample size
a<-array(0,m)
ll<-array(0,m)
ul<-array(0,m)
i<-0
while (i<m) {i<-i+1
a[i]<-mean(rnorm(n))
ll[i]<-a[i]+qnorm(0.025)*sqrt(1/n)
ul[i]<-a[i]+qnorm(0.975)*sqrt(1/n)}
interval.plot(ll,ul)
#
# Example 8.2
xbar <- 77/30
Z <- qnorm(.975)
sigma <- 2.5
n <- 30
# 95% CI
round(c(xbar-Z*sigma/sqrt(n),xbar+Z*sigma/sqrt(n)),2)
# Or Using zsum.test
round(zsum.test(mean.x=xbar,sigma.x=sigma,n.x=n)$conf,2)

#
# Example 8.3
attach(Grocery)
qqnorm(groceries)
qqline(groceries)
#
mean(groceries)
qnorm(0.03/2)
round(c(mean(groceries)-qnorm(1-0.03/2)*sqrt(900/30),
        mean(groceries)+qnorm(1-0.03/2)*sqrt(900/30)),2)
# OR
xbar <- mean(groceries)
Z <- qnorm(.985)
sigma <- 30
n <- 30
# 97% CI
round(c(xbar-Z*sigma/sqrt(n),xbar+Z*sigma/sqrt(n)),2)
# OR using z.test() (PASWR)
round(z.test(x=groceries,sigma.x=30,conf.level=.97)$conf,2)
# part b.
z.test(x=groceries,sigma.x=30,conf.level=.99,alternative="greater")$conf
detach(Grocery)
#
# Example 8.4
n <- ceiling(qnorm(.975)^2*.1^2/.02^2)
n
# OR using nsize() (PASWR)
nsize(b=0.02, sigma=0.1, conf.level=.95, type="mu")
#
# Example 8.5
n <- ceiling(qnorm(.975)^2*5^2)
n
# Example 8.6
# a.
f <- 0 # f= number of failures
p <- 0 # p = probability
s <- 100 # s = number of successes
while(p < 0.95)
{
f <- f + 1
p <- pnbinom(f, s, 0.03)
}
ans <- c(f + s, p) # f+s= Containers
names(ans) <- c("Containers", "Probability")
ans
# b.
f <- function(x)
{
0.0009*x^2 - 6.07873131*x +10000
}
curve(f,2800,4000)
abline(h=0,col="red",lty=2)
uniroot(f,c(3000,4000))
# c.
n <- 0 # Number of containers
p <- 0 # Probability
while(p < 0.95)
{
n <- n + 1
p <- 1 - pbinom(99, n, 0.03)}
ans <- c(n, p)
names(ans) <- c("Containers", "Probability")
ans
#
#
# Example 8.7
xbar <- 61.9/12
S <- sqrt((450 - 12*xbar^2)/11)
ct <- qt(0.95,11)
round(c(xbar-ct*S/sqrt(12), xbar+ct*S/sqrt(12)),2)
# OR
round(tsum.test(mean.x=xbar,s.x=S,n.x=12,conf.level=.90)$conf,2)
#
#
# Example 8.8
attach(House)
MEAN<-mean(Price)
CT<-qt(.975,13)
ST<-sd(Price) #stdev(Price) in S-PLUS
round(c(MEAN-CT*ST/sqrt(14), MEAN+CT*ST/sqrt(14)),2)
# OR
round(t.test(Price)$conf,2)
detach(House)
#
#
# Example 8.9
round(qnorm(0.975), 2)
round(c((4-4.41) + qnorm(0.025)*3*sqrt(1/15 + 1/22),
        (4-4.41) + qnorm(0.975)*3*sqrt(1/15 + 1/22)), 2)
# Or
xbar <- 60/15
n.x <- 15
sigma <- 3
ybar <- 97/22
n.y <- 22
Z <- qnorm(.975)
round(c((xbar-ybar) - Z*sigma*sqrt(1/n.x + 1/n.y),
(xbar-ybar) + Z*sigma*sqrt(1/n.x + 1/n.y)),2)
# Or
round(zsum.test(mean.x=xbar, mean.y=ybar, sigma.x=sigma,
sigma.y=sigma, n.x=15, n.y=22)$conf, 2)
#
#
# Example 8.10
# Figure 8.5
attach(Apple)
par(pty = "s")
Altblue <- "#A9E2FF"
Adkblue <- "#0080FF"
fresh <- qqnorm(Fresh)
old <- qqnorm(Warehouse)
plot(fresh,type="n",ylab="Sample Quantiles",
xlab="Theoretical Quantiles")
qqline(Fresh, col = Altblue)
qqline(Warehouse, col = Adkblue)
points(fresh, col = Altblue, pch = 16, cex = 1.2)
points(old, col = Adkblue, pch = 17)
legend(-1.75, 9.45, c("Fresh", "Warehouse"),
col = c(Altblue, Adkblue), text.col=c("black","black"),
pch=c(16,17), lty=c(1,1), bg="gray95", cex=0.75)
title("Q-Q Normal Plots")
#
mean.fresh <- mean(Fresh)
mean.fresh
mean.old <- mean(Warehouse)
mean.old
round(c(mean.fresh - mean.old - qnorm(0.975)*1.5*sqrt(2/17),
        mean.fresh - mean.old + qnorm(0.975)*1.5*sqrt(2/17)),2)
# 95% CI using z.test()
round(z.test(x=Fresh, y=Warehouse, sigma.x=1.5,
sigma.y=1.5)$conf,2)
detach(Apple)
#
#
# Example 8.11
#
round(c((8.4 - 8.8) - qnorm(0.985)*sqrt((4.5)^2/50 + (6)^2/46),
        (8.4 - 8.8) + qnorm(0.985)*sqrt((4.5)^2/50 + (6)^2/46)),2)
# OR without the rounding in ybar using zsum.test()
round(zsum.test(mean.x=420/50, mean.y=405/46, sigma.x= 4.5,
sigma.y=6, n.x=50, n.y=46, conf.level=0.97)$conf,2)
#
#
# Example 8.12
attach(Calculus)
str(Calculus) # str ONLY works in R
mean(Yes.Calculus)
mean(No.Calculus)
round(c(mean(Yes.Calculus) - mean(No.Calculus)
- qnorm(0.975)*sqrt(25/18+144/18),
mean(Yes.Calculus) - mean(No.Calculus)
+ qnorm(0.975)*sqrt(25/18+144/18)),2)
# OR using z.test()
round(z.test(x=Yes.Calculus, y=No.Calculus,
sigma.x=5, sigma.y=12)$conf, 2)
detach(Calculus)
#
#
# Example 8.13
#
xbar <- 53/15
sx <- sqrt((222-15*xbar^2)/14)
nx <- 15
ybar <- 77/11
sy <- sqrt((560-11*ybar^2)/10)
ny <- 11
sp <- sqrt(((nx-1)*sx^2 +(ny-1)*sy^2)/(nx+ny-2))
sp
round(c((xbar-ybar)-qt(0.975,24)*sp*sqrt(1/nx+1/ny),
        (xbar-ybar)+qt(0.975,24)*sp*sqrt(1/nx+1/ny)),2)
# OR using tsum.test()
round(tsum.test(mean.x=xbar, mean.y=ybar, s.x=sx, s.y=sy,
n.x=nx, n.y=ny, var.equal=TRUE)$conf, 2)
#
#
# Example 8.14
attach(Apple)
MX <- mean(Fresh)
MX
VX <- var(Fresh)
VX
MY <- mean(Warehouse)
MY
VY <- var(Warehouse)
VY
SP <- sqrt((16*VX + 16*VY)/32)
SP
CT <- qt(.975,32)
CT
c( (MX-MY) - CT*SP*sqrt(1/17 + 1/17) ,
   (MX-MY) + CT*SP*sqrt(1/17 + 1/17) )
t.test(Fresh, Warehouse, var.equal=TRUE)$conf
detach(Apple)
#
#
# Example 8.15
xbar <- 63/15
xbar
nx <- 15
sx <- sqrt((338-15*xbar^2)/14)
sx
ybar <- 66.4/11
ybar
ny <- 11
sy <- sqrt((486-11*ybar^2)/10)
sy
nu <- (sx^2/15+sy^2/11)^2 / (((sx^2/15)^2/14) + ((sy^2/11)^2/10))
nu
# Note the book uses 18 degrees of freedom instead of 18.38828
round(c((xbar-ybar) - qt(.975,18)*sqrt(sx^2/15+sy^2/11),
        (xbar-ybar) + qt(.975,18)*sqrt(sx^2/15+sy^2/11) ), 2)
# Or using nu degrees of freedom
round(c((xbar-ybar) - qt(.975,nu)*sqrt(sx^2/15+sy^2/11),
        (xbar-ybar) + qt(.975,nu)*sqrt(sx^2/15+sy^2/11) ), 2)
# Or
round(tsum.test(mean.x=xbar, mean.y=ybar, s.x=sx,
s.y=sy, n.x=nx, n.y=ny, var.equal=FALSE)$conf, 2)
#
#
# Example 8.16
attach(Calculus)
t.test(Yes.Calculus, No.Calculus, var.equal=FALSE)$conf
detach(Calculus)
#
#
# Example 8.17
#
attach(Sundig)
# Figure 8.7
round(qt(0.025,4), 2)
qqnorm(d)
qqline(d)
# Similar to Figure 8.8
ntester(d)
#
t.test(SUN, DIGITAL, paired=TRUE)$conf
detach(Sundig)
#
#
# Example 8.18
#
# Uses Example 8.15 information:
xbar <- 63/15
nx <- 15
sx <- sqrt((338-15*xbar^2)/14)
sx^2
ybar <- 66.4/11
ny <- 11
sy <- sqrt((486-11*ybar^2)/10)
sy^2
# 80% CI
round(c((nx-1)*sx^2/qchisq(.9,nx-1),
        (nx-1)*sx^2/qchisq(.1,nx-1)), 2)
#
qchisq(0.1,14)
qchisq(0.9,14)
round(c(14*5.24/qchisq(0.9,14), 14*5.24/qchisq(0.1,14)), 2)
#
# Example 8.19
# a.
library(lattice)
attach(barley)
n <- length(yield[year==1932])
n
mean(yield[year==1932])
var(yield[year==1932])
sd(yield[year==1932]) #S-PLUS: stdev(yield[year==1932]
qt(.975,n-1)
#
# To construct the confidence interval directly with S
t.test(yield[year==1932])$conf
#
# b.
s2<- var(yield[year==1932])
s2
ChiL <- qchisq(.025,59)
ChiL
ChiU <- qchisq(.975,59)
ChiU
n <- length(yield[year==1932])
n
round(c((n-1)*s2/ChiU, (n-1)*s2/ChiL),2)
detach(barley)
#
#
# Example 8.20
# Note: this code does not round the
# numbers until the very end...hence the slight difference
# between what appears in the book.
# information from 8.13:
xbar <- 53/15
sx <- sqrt((222-15*xbar^2)/14)
nx <- 15
ybar <- 77/11
sy <- sqrt((560-11*ybar^2)/10)
ny <- 11
# 90% CI
round( c( qf(.05,nx-1,ny-1)*sx^2/sy^2,
          qf(.95,nx-1,ny-1)*sx^2/sy^2),2)
#
#
# Example 8.21
attach(Apple)
var(Warehouse)
var(Fresh)
round(var(Warehouse)/var(Fresh),2)
round(qf(0.025,16,16), 2)
round(qf(0.975,16,16), 2)
var.test(Warehouse, Fresh)$conf
detach(Apple)
#
#
# Example 8.22
round(c(2-1.64*sqrt(0.02), 2+1.64*sqrt(0.02)), 2)
#
#
# Example 8.23
#
# a.
p <- 26/40
p
z <- qnorm(.975)
z
n <- 40
round(z,2) # z_(0.975)
round(c(p - z*sqrt(p*(1 - p)/n),
        p + z*sqrt(p*(1 - p)/n)),3)
# b.
round(prop.test(26,40,correct=FALSE)$conf, 3)
# c.
round(prop.test(26,40,correct=TRUE)$conf, 3)
#
#
# Example 8.24
round(prop.test(157, 200, conf.level=0.90, correct=FALSE)$conf,2)
round(prop.test(157, 200, conf.level=0.95, correct=FALSE)$conf,2)
round(prop.test(157, 200, conf.level=0.99, correct=FALSE)$conf,2)
#
#
# Example 8.25
# a.
p <- 400/2000
n <- 2000
z <- qnorm(.975)
round(z,2)
round(c(p-z*sqrt(p*(1-p)/n), p+z*sqrt(p*(1-p)/n)), 3)
# b.
p <- seq(0, 1, 0.001)
plot(p, p*(1-p), type="l")
n <- ceiling( (qnorm(.975)/0.015)^2*.5*.5 )
n
# OR
nsize(b=0.015,p=0.5,type="pi")
#
#
# Example 8.26
round(prop.test(c(70,90),c(1400,2000), correct=FALSE)$conf,4)
round(prop.test(c(70,90),c(1400,2000), correct=TRUE)$conf,4)
#
#
# Example 8.27
attach(Soccer)
M <- mean(Goals,na.rm=T)
M
z <- qnorm(.95)
z
round(c(M-z*sqrt(M/232),M+z*sqrt(M/232)),2)
detach(Soccer)